<p>Using Nextflow to submit jobs this way makes little sense, since:</p>
<blockquote>
<p><a href="https://slurm.schedmd.com/sbatch.html" rel="nofollow noreferrer">sbatch</a> exits immediately after the script is successfully transferred
to the Slurm controller and assigned a Slurm job ID. The batch script
is not necessarily granted resources immediately, it may sit in the
queue of pending jobs for some time before its required resources
become available.</p>
</blockquote>
<p>Assuming your job script is basically a shell script, make it executable (e.g. <code>chmod +x PRE_PHASING_QC.JOB</code>) and move it into a folder called <code>bin</code> in the root directory of your project repository. <a href="https://www.nextflow.io/docs/latest/sharing.html#third-party-scripts" rel="nofollow noreferrer">Nextflow automatically adds this folder to your $PATH</a>, making scripts in this folder accessible in your pipeline without needing localize them and without needing to specify an absolute path to invoke them. Your code then becomes:</p>
<pre><code>process prephasing {

    &quot;&quot;&quot;
    PRE_PHASING_QC.JOB [options] &lt;args&gt;
    &quot;&quot;&quot;
}
</code></pre>
<p>Finally, add <code>process.executor = 'slurm'</code> to your <code>nextflow.config</code> file. You may need to configure some other options, like resource requests, but that should get your script integrated into your pipeline without much refactoring, if any.</p>
<hr />
<hr />
<p>Contents of <code>PRE_PHASING_QC.sh</code>:</p>
<pre><code>#!/bin/bash

set -euo pipefail

hrc_sites_file=&quot;<span class="math-container">${1}"
input_prefix="$</span>{2}&quot;
chrom=&quot;chr${3}&quot;

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs)

plink --bfile &quot;<span class="math-container">${input_prefix}" --chr "$</span>{3}&quot; --out &quot;<span class="math-container">${chrom}" --make-bed
plink --bfile "$</span>{chrom}&quot; --list-duplicate-vars ids-only suppress first --out &quot;<span class="math-container">${chrom}.step1a"
plink --bfile "$</span>{chrom}&quot; --exclude &quot;<span class="math-container">${chrom}.step1a.dupvar" --make-bed --out "$</span>{chrom}.step1a&quot;
plink --bfile &quot;<span class="math-container">${chrom}.step1a" --snps-only just-acgt --make-bed --out "$</span>{chrom}.step1b&quot;

awkcmd='
(<span class="math-container">$5 == "T" &amp;&amp; $</span>6 == &quot;A&quot;) || (<span class="math-container">$5 == "A" &amp;&amp; $</span>6 == &quot;T&quot;) ||
(<span class="math-container">$5 == "C" &amp;&amp; $</span>6 == &quot;G&quot;) || (<span class="math-container">$5 == "G" &amp;&amp; $</span>6 == &quot;C&quot;)
{
    print <span class="math-container">$2, "ambig"
    next
}
{
    print $</span>2
}'

awk &quot;<span class="math-container">${awkcmd}" "$</span>{chrom}.step1b&quot; | grep &quot;ambig&quot; &gt; &quot;<span class="math-container">${chrom}.step1b.snplist.txt"
plink --bfile "$</span>{chrom}.step1b&quot; --exclude &quot;<span class="math-container">${chrom}.step1b.snplist.txt" --make-bed --out "$</span>{chrom}.step1c&quot;

# Step 2. Zero out Mendelian errors
plink --bfile &quot;<span class="math-container">${chrom}.step1c" --me 1 1 --set-me-missing --make-bed --out "$</span>{chrom}.step2&quot;

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk.
plink --bfile &quot;<span class="math-container">${chrom}.step2" --hwe 1e-6 --geno 0.02 --make-bed --out "$</span>{chrom}.step3&quot;
plink --bfile &quot;<span class="math-container">${chrom}.step3" --freq --make-bed --out "$</span>{chrom}.step4&quot;

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s)
HRC-1000G-check-bim.pl -b &quot;<span class="math-container">${chrom}.step4.bim" -f "$</span>{chrom}.step4.frq&quot; -r &quot;<span class="math-container">${hrc_sites_file}" -h
plink --bfile "$</span>{chrom}.step4&quot; --exclude &quot;Exclude-<span class="math-container">${chrom}.step4-HRC.txt" --make-bed --out "$</span>{chrom}.step5&quot;
plink --bfile &quot;<span class="math-container">${chrom}.step5" --update-chr "Chromosome-$</span>{chrom}.step4-HRC.txt&quot; --make-bed --out &quot;<span class="math-container">${chrom}.step6"
plink --bfile "$</span>{chrom}.step6&quot; --update-map &quot;Position-<span class="math-container">${chrom}.step4-HRC.txt" --make-bed --out "$</span>{chrom}.step7&quot;
plink --bfile &quot;<span class="math-container">${chrom}.step7" --flip "Strand-Flip-$</span>{chrom}.step4-HRC.txt&quot; --make-bed --out &quot;<span class="math-container">${chrom}.step8"
plink --bfile "$</span>{chrom}.step8&quot; --a2-allele &quot;Force-Allele1-<span class="math-container">${chrom}.step4-HRC.txt" --make-bed --out "$</span>{chrom}.step9&quot;
plink --bfile &quot;<span class="math-container">${chrom}.step9" --update-name "ID-$</span>{chrom}.step4-HRC.txt&quot; --real-ref-alleles --make-bed --out &quot;${chrom}.step10&quot;

</code></pre>
<p>Make <code>PRE_PHASING_QC.sh</code> executable (<code>chmod +x PRE_PHASING_QC.sh</code>), and replace your job script with the following. This file won't be needed by your Nextflow process, but shows you how the above shell script can be invoked:</p>
<pre><code>#!/bin/bash
#SBATCH --job-name=prephase
#SBATCH --account=p697_norment
#SBATCH --time=1:00:00
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=7600M
#SBATCH --cpus-per-task=7
#SBATCH --array=1-22

set -e

. /cluster/bin/jobsetup

export PATH=&quot;/cluster/projects/p697/github/norment/moba_qc_imputation/software:${PATH}&quot;

module load plink

cd /path/to/hce-ec-eur-fin-pass-qc

/path/to/PRE_PHASING_QC.sh \
    /cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz \
    hce-ec-eur-fin-pass-qc \
    &quot;${SLURM_ARRAY_TASK_ID}&quot;
</code></pre>
<p>Copy <code>PRE_PHASING_QC.sh</code> and the Oxford Perl script, <code>HRC-1000G-check-bim.pl</code>, into your project's <code>bin/</code> directory. You may also need to make the Perl script executable. Then your Nextflow process becomes:</p>
<pre><code>hrc_sites_file = file(params.hrc_sites_file)

Channel
    .fromFilePairs( &quot;${params.input_dir}/*.{bed,bim,fam}&quot;, size:3 )
    .set { plink_data }

Channel
    .of(1..22)
    .set { chroms }


process prephasing {

    input:
    tuple val(prefix), path(bfiles) from plink_data
    each chrom from chroms
    path hrc_sites_file

    &quot;&quot;&quot;
    PRE_PHASING_QC.sh &quot;<span class="math-container">${hrc_sites_file}" "$</span>{prefix}&quot; &quot;${chrom}&quot;
    &quot;&quot;&quot;
}
</code></pre>
<p>Once the above runs successfully, you can declare your process outputs. Please note that none of the above was tested; hopefully not too many typos, but YMMV.</p>

CC BY-SA 4.0
2021-03-28T22:59:56.653Z
Thanks `@Steve` . I think you understand my explanation and your answer makes sense. I have updated the question with new code, error and expected output. Hope it helps to debug
zillur rahman
CC BY-SA 4.0
2021-03-28T23:48:32.140Z
@zillurrahman: You get the error `Missing value declared as output parameter: chr` because you've defined 'chr' in your output block, but this value wasn't defined in your input block (or globally in your script block). You didn't show what your input block looks like. I suspect it should look something like: `tuple val(chr), val(num), val(steps), path(bfiles)`. I don't know what you job script does, but if your job script references absolute paths to files, you'll need to refactor it so that they can be passed in as variables on the command line.
Steve
CC BY-SA 4.0
2021-03-29T07:45:19.477Z
My job script takes bfiles as input (e.g `plink --bfiles input_bfiles --other argument --output`). I just tried to modify this as : `plink --bfiles $1` and tried to use input channel from nextflow. `""" 
PRE_PHASING_QC.job "${bfiles}"
""""` In this case nothing happens, no error or success.
zillur rahman
CC BY-SA 4.0
2021-03-29T10:02:51.523Z
From our previous discussions, I thought plink took a prefix for it's '--bfile' option? Note this isn't the same as '--bfiles'. If I'm correct, I think your first task will be to work out how to supply a prefix on the command line so that it gets passed down to your plink command. In your NF process, I think you'll want: `PRE_PHASING_QC.job "${prefix}"`. And in your script I think you'll want: `plink --bfile "$1"`
Steve
CC BY-SA 4.0
2021-03-30T12:31:03.447Z
Thanks `@Steve` . I tried `nextflow run prephasing1.nf --hrc_sites_file "path/to/hrc_sites.tab.gz" "path/to/hce-ec-eur-fin-pass-qc" "1"` It give me `Missing 'bind' declaration in input parameter` . Should we use `tuple(val),val(input_prefix),path(bfiles) from some_channel` ? How could I declare input channel in this particular case? It seems I could run chromosome by chromosome? That if I put "1" in nf run parameter, it will only run for chromosome 1? Should we declare it in `in_prefix=file(params.input_prefix)` and `params.chr=1` then we could change/call the arguments in command line?
zillur rahman
CC BY-SA 4.0
2021-03-30T12:50:25.573Z
Thanks `@Steve` So far I declared input_prefix and chrom as hrc_sites_file was declared in your script and specify them in command line. I got output for chromosome 1, chr1.step1a, chr1.step1b but in the next step it shows that no such file as chr1.step1b but it exists in my work/ directory. Any thought? When I use chrom 2, I got chr2.step1a, chr2.step1b and same error in the next step.
zillur rahman
CC BY-SA 4.0
2021-03-30T14:19:29.910Z
No worries @zillurrahman. I just updated the NF code to process all chromosomes like what we have done previously. I think you'd rather avoid running each chromosome one at a time. In your first comment above, you get the `Missing 'bind' declaration in input parameter` error since you're missing parameters for the `"path/to/hce-ec-eur-fin-pass-qc"` and `"1"` options.
Steve
CC BY-SA 4.0
2021-03-30T20:40:25.327Z
Ohh! There was some typo in the PRE_PHASING_QC.sh script. I ran without the slurm, identfied those typos, edited them and now everything is Ok. Now, How could I declare the output files so that I get all the outputs in a single directory? So far, I tried `tuple val("${chrom}"), val(steps), path("${chrom}.${steps}.{bed,bim,fam}") into qc1_bfile` and `tuple val(prefix), val("${chrom}"), path("${prefix}-${chrom}.step4-HRC.txt") into qc1_hrc` But got error:  `Missing value decalred as output parameter: steps` How could I decalre steps so that all outputs from each step gathered in a folder?
zillur rahman
CC BY-SA 4.0
2021-03-30T21:49:10.587Z
Got it! I did as following: `tuple val("${chrom}"), path("chr${chrom}*.{bed,bim,fam}") into qc1_chr` and it works. Only thing is, file transfer taking too long. After finishing the process `Waiting files transfer to complete (2 files)` it is taking too much time. If I check output folder, I can see all my desired outputs but the waiting files transfer to complete message is still on the screen and I can not do anything in this tab. Any suggestion? However, Thanks for helping all the way!
zillur rahman
CC BY-SA 4.0
2021-03-31T03:33:15.413Z
Well done getting it working @zillurrahman. Is the issue that the file copying is taking too long, or is the issue just that while you're waiting you cannot use the current terminal? If it's the latter, it's not really a problem since the job can be run non-interactively (e.g. submitted as a job to your Slurm cluster), in the background (e.g. using screen or nohup) or can be backgrounded after i'ts started (e.g. using <ctrl-z>, then typing 'bg') to background the process.
Steve
CC BY-SA 4.0
2021-03-31T08:42:00.657Z
Thanks `@Steve`. I have managed to run it using `sbatch`. It is ok now.
zillur rahman
CC BY-SA 4.0
2021-04-02T23:07:57.267Z
Hi, This process creates multiple output files `(chr1-22.step1-10.bed|bim|fam)` which are saved into qc1_chr  output channel. Now, say I need only those files which has string `step10`. How could I get that subset into a new channel? I need it into a new channel because I need to merge them with elements of other channels. Did my explanation make sense?
zillur rahman
CC BY-SA 4.0
2021-04-03T11:16:53.247Z
@zillurrahman you could use separate output channels for each step, but you might not need the other step outputs. Another way would be to select just the files you need from the qc1_chr output channel using the [filter](https://www.nextflow.io/docs/latest/operator.html#filter) operator. There's lots of ways to do this - for example: `qc1_chr.filter { bfile -> ["bed", "bim", "fam"].any { bfile.name.endsWith("step10.${it}") } }.view()`. HTH.
Steve
CC BY-SA 4.0
2021-04-03T16:24:35.637Z
Thanks `@Steve`. It gives me `No such variable: Exception evaluating property 'name' for java.util.ArrayList, Reason: groovy.lang.MissingPropertyException` . If I remove `name` from the above snippet, I got `Unknown method invocation 'endsWith' on ArrayList type`. Any idea?
zillur rahman
CC BY-SA 4.0
2021-04-04T10:58:44.317Z
Sorry @zillurrahman I didn't see your qc1_chr declaration above. I think you want: `qc1_chr.filter { name, bfiles -> bfiles.collect { it.getBaseName() }.every { it.endsWith("step10") } }` HTH.
Steve
CC BY-SA 4.0
2021-04-04T17:48:27.240Z
Thanks `@Steve`. It seems ok now, getting no error. Another issue. I merged other channels with it (each channel for 22 files, 1 for each chromosome) `shape_para=shape_in.merge(shape_map).merge(shape_ref).merge(shape_legend).merge(shape_sample)
process phasing {
input:
each file(penta) from shape_para
output:
path("*step10*") into phased
script:
def (in, map, ref, legend, sample) = penta
"""
shapeit2 -T 32 --input-bed "${in}" --input-map "${map}" --input-ref "${ref}" "${legend}" "${sample}" --duohmm -O "${in}.phased" --seed 54321
"""
}` I am getting no output or error!
zillur rahman
CC BY-SA 4.0
2021-04-05T10:14:54.583Z
I'm not sure exactly @zillurrahman, sorry. Your input declaration doesn't look right, but I'd still expect some output or log on stderr if your process is executed. If you can [ask a new question](https://bioinformatics.stackexchange.com/questions/ask), I'll see if I can add an answer. Please also show the contents of `shape_in`, `shape_ref`, `shape_legend` and `shape_sample`. Cheers.
Steve
CC BY-SA 4.0
2021-04-05T10:25:18.860Z
Thanks `@Steve` I posted a new question with code for each channel here; https://bioinformatics.stackexchange.com/questions/15680/merge-multiple-channels-each-emits-multiple-files-in-nextflow
zillur rahman
CC BY-SA 4.0
2021-03-28T12:28:23.003Z
2021-03-30T14:15:10.390Z
2021-03-30T14:15:10.390Z
Steve
<p><a href="https://i.sstatic.net/2dD9s.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/2dD9s.png" alt="enter image description here" /></a><a href="https://i.sstatic.net/A23k1.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/A23k1.jpg" alt="enter image description here" /></a>I was trying to run a slurm job script inside nextflow process. Is it possible? For instance, I have a job script 'PRE_PHASING_QC.JOB' and outside nextflow it runs perfectly, I want to integrate it in a nextflow script. May be my example code could be an explanation of what I want.</p>
<pre><code>prephase='PRE_PHASING_QC.JOB'
process prephasing {
input:
file(prephasing) from prephase
&quot;&quot;&quot;
sbatch &quot;${prephasing}
&quot;&quot;&quot;
}
</code></pre>
<p>The error:</p>
<pre><code>sbatch: error: This does not look like a batch script
...
</code></pre>
<p>I tried to make my job script executable, put it in bin and then include it in the nextflow process according to @Steve. So, the script runs without any error but got no output. The original job script takes 3 files (bed,bim,fam) as input and produce multiple outputs. I guess that I have to declare outputs inside nextflow process. So far I tried as follows:</p>
<pre><code>process prephasing {
publishDir params.out, mode: 'copy'
output:
tuple val(chr),val(num),val(steps),path(&quot;<span class="math-container">${chr}$</span>{num}.${steps}.{bed,bim,fam,log}&quot;) into out1
script:
&quot;&quot;&quot;
PRE_PHASING_QC.job
&quot;&quot;&quot;
</code></pre>
<p>It gives:</p>
<pre><code>Error executing process &gt; 'prephasing'
Caused by:
Missing value declared as output parameter: chr
</code></pre>
<p>My expected output would be chr{1-22}.step{1-10}.{bed,bim,fam,log}</p>
<pre><code>#!/bin/bash 
#SBATCH --job-name=prephase 
#SBATCH --account=p697_norment 
#SBATCH --time=1:00:00 
#SBATCH --ntasks=1 
#SBATCH --mem-per-cpu=7600M 
#SBATCH --cpus-per-task=7 
#SBATCH --array=1-22 

source /cluster/bin/jobsetup 

set -o errexit 

module load plink 

export prefix=chr${SLURM_ARRAY_TASK_ID} 
export hrc_sites_file=/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz 
export HRC_1000G_check_bim=/cluster/projects/p697/github/norment/moba_qc_imputation/software/HRC-1000G-check-bim.pl

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs) 

plink --bfile hce-ec-eur-fin-pass-qc --chr <span class="math-container">${SLURM_ARRAY_TASK_ID} --out $</span>prefix --make-bed || true
plink --bfile <span class="math-container">$prefix --list-duplicate-vars ids-only suppress first --out $</span>prefix.stepla || true 
plink --bfile <span class="math-container">$prefix --exclude $</span>prefix.stepla.dupvar --make-bed --out <span class="math-container">$prefix.stepla || true 
plink --bfile $</span>prefix.stepla --snps-only just-acgt --make-bed --out $prefix.steplb || true 

awk '{ if ((<span class="math-container">$5=="T" &amp;&amp; $</span>6==&quot;A&quot;)||(<span class="math-container">$5=="A" &amp;&amp; $</span>6==&quot;T&quot;)ii(<span class="math-container">$5=="C" &amp;&amp; $</span>6==&quot;G&quot;)||(<span class="math-container">$5=="G" &amp;&amp; $</span>6==&quot;C&quot;)) print <span class="math-container">$2, "ambig" ; else print $</span>2 ;}' <span class="math-container">$prefix. steplb.bim | grep ambig &gt; $</span>prefix.steplb.snplist.txt || true 
plink --bfile <span class="math-container">$prefix.steplb --exclude $</span>prefix.steplb.snplist.txt --make-bed --out $prefix.steplc || true 

# Step 2. Zero out Mendelian errors
plink --bfile <span class="math-container">$prefix.steplc --me 1 1 --set-me-missing --make-bed --out $</span>prefix.step2 || true 

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies 
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk. 
plink --bfile <span class="math-container">$prefix.step2 --hwe le-6 --geno 0.02 --make-bed --out $</span>prefix.step3 || true 
plink --bfile <span class="math-container">$prefix.step3 --freq --make-bed --out $</span>prefix.step4 || true 

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s) 
perl <span class="math-container">${HRC-IRC_1000G_check_bim} -b $</span>prefix.step4.bim -f <span class="math-container">$prefix.step4.frq -r $</span>{hrc_sites_file} -h || true 
plink --bfile <span class="math-container">$prefix.step4 --exclude Exclude-$</span>prefix.step4-HRC.txt --make-bed --out <span class="math-container">$prefix.step5 || true 
plink --bfile $</span>prefix.step5 --update-chr Chromosome-<span class="math-container">$prefix.step4-HRC.txt --make-bed --out $</span>prefix.step6 || true 
plink --bfile <span class="math-container">$prefix.step6 --update-map Position-$</span>prefix.step4-HRC.txt --make-bed --out <span class="math-container">$prefix.step7 || true 
plink --bfile $</span>prefix.step7 --flip Strand-Flip-<span class="math-container">$prefix.step4-HRC.txt --make-bed --out $</span>prefix.step8 || true 
plink --bfile <span class="math-container">$prefix.step8 --a2-allele Force-Allele1-$</span>prefix.step4-HRC.txt --make-bed --out <span class="math-container">$prefix.step9 || true 
plink --bfile $</span>prefix.step9 --update-name ID-<span class="math-container">$prefix.step4-HRC.txt --real-ref-alleles --make-bed --out $</span>prefix.step10 || true
</code></pre>

CC BY-SA 4.0
2021-03-28T23:59:37.780Z
When you say your script ran without error but you got no output files, are you looking in the workDir or publishDir? When I'm developing, I usually only write the output block declaration after I've tested the tool/script a couple of times. This is because your task process can still run successfully without an output block declaration. I can then look in the workdir to see what files have been created, etc, then finally declare them in the output block. Then I start on the next (downstream) process(es) and repeat.
Steve
CC BY-SA 4.0
2021-03-29T07:37:24.377Z
The publishDir is empty and in the workDir there are only log files (chr.log, chr.step1.log, chr.step2.log etc.). It is strange!
zillur rahman
CC BY-SA 4.0
2021-03-29T10:08:42.363Z
So only files declared in your output block will get published into your publishDir. Ignore output in the publishDir while testing and developing. Or better, comment out the publishDir directive until you get your script block working in some capacity. In any case, your job script will show you what's going on. Happy to take a look if you need help.
Steve
CC BY-SA 4.0
2021-03-29T13:15:20.317Z
Thank you `@Steve`. I added the job file as a png (not sure how to add text file). In line 21 `hce-ec-eur-fin-pass-qc` is the only input which create chr(1-22).(bed,bim,fam) files (for 22 chromosomes) which are used in the subsequent commands. I tried to replace `hce-ec-eur-fin-pass-qc` with `$1` but got no output. Maybe you could give a better idea.
zillur rahman
CC BY-SA 4.0
2021-03-29T22:45:33.393Z
Please correct me if I'm wrong, but it looks like this script kicks off a job array. This isn't a good fit inclusion "as-is" into a Nextflow workflow, since it won't be portable at all. I think breaking it up into two scripts is the best option here. I'll see if I can give an example of this. Also, you can paste text in a code-block. I've updated your question with the text. Click [edit](https://bioinformatics.stackexchange.com/posts/15633/edit) to see the changes I've made.
Steve
CC BY-SA 4.0
2021-03-29T22:57:11.157Z
See also: https://bioinformatics.stackexchange.com/editing-help
Steve
CC BY-SA 4.0
2021-03-30T03:12:51.463Z
Just realized my edit hasn't been approved yet, but you should be able to see it.
Steve
CC BY-SA 4.0
2021-03-30T06:49:39.533Z
Thanks `@Steve` . Yes, it starts a job array 1-22 for each chromosome. The initial input is one single bfile (bed,bim,fam), in step1a it creates 22 bfiles (bed,bim,fam) for each chromosome. The rest of steps operate on those 22 files to create final output. Any relatable example to break the job into two would be great.
zillur rahman
CC BY-SA 4.0
2021-03-30T10:32:07.867Z
No worries at all @zillurrahman. I've updated my answer with an example on how to do this. I also refactored some of your code, hopefully it still works...
Steve
<nextflow>
Nextflow: run sbatch (slurm job script) inside a process
# Nextflow: run sbatch (slurm job script) inside a process

(Asked by: username_0 on 2021-03-28)

I was trying to run a slurm job script inside nextflow process. Is it possible? For instance, I have a job script 'PRE_PHASING_QC.JOB' and outside nextflow it runs perfectly, I want to integrate it in a nextflow script. May be my example code could be an explanation of what I want.
prephase='PRE_PHASING_QC.JOB'
process prephasing {
input:
file(prephasing) from prephase
"""
sbatch "${prephasing}
"""
}

The error:
sbatch: error: This does not look like a batch script
...

I tried to make my job script executable, put it in bin and then include it in the nextflow process according to @Steve. So, the script runs without any error but got no output. The original job script takes 3 files (bed,bim,fam) as input and produce multiple outputs. I guess that I have to declare outputs inside nextflow process. So far I tried as follows:
process prephasing {
publishDir params.out, mode: 'copy'
output:
tuple val(chr),val(num),val(steps),path("${chr}${num}.${steps}.{bed,bim,fam,log}") into out1
script:
"""
PRE_PHASING_QC.job
"""

It gives:
Error executing process > 'prephasing'
Caused by:
Missing value declared as output parameter: chr

My expected output would be chr{1-22}.step{1-10}.{bed,bim,fam,log}
#!/bin/bash 
#SBATCH --job-name=prephase 
#SBATCH --account=p697_norment 
#SBATCH --time=1:00:00 
#SBATCH --ntasks=1 
#SBATCH --mem-per-cpu=7600M 
#SBATCH --cpus-per-task=7 
#SBATCH --array=1-22 

source /cluster/bin/jobsetup 

set -o errexit 

module load plink 

export prefix=chr${SLURM_ARRAY_TASK_ID} 
export hrc_sites_file=/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz 
export HRC_1000G_check_bim=/cluster/projects/p697/github/norment/moba_qc_imputation/software/HRC-1000G-check-bim.pl

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs) 

plink --bfile hce-ec-eur-fin-pass-qc --chr ${SLURM_ARRAY_TASK_ID} --out $prefix --make-bed || true
plink --bfile $prefix --list-duplicate-vars ids-only suppress first --out $prefix.stepla || true 
plink --bfile $prefix --exclude $prefix.stepla.dupvar --make-bed --out $prefix.stepla || true 
plink --bfile $prefix.stepla --snps-only just-acgt --make-bed --out $prefix.steplb || true 

awk '{ if (($5=="T" && $6=="A")||($5=="A" && $6=="T")ii($5=="C" && $6=="G")||($5=="G" && $6=="C")) print $2, "ambig" ; else print $2 ;}' $prefix. steplb.bim | grep ambig > $prefix.steplb.snplist.txt || true 
plink --bfile $prefix.steplb --exclude $prefix.steplb.snplist.txt --make-bed --out $prefix.steplc || true 

# Step 2. Zero out Mendelian errors
plink --bfile $prefix.steplc --me 1 1 --set-me-missing --make-bed --out $prefix.step2 || true 

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies 
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk. 
plink --bfile $prefix.step2 --hwe le-6 --geno 0.02 --make-bed --out $prefix.step3 || true 
plink --bfile $prefix.step3 --freq --make-bed --out $prefix.step4 || true 

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s) 
perl ${HRC-IRC_1000G_check_bim} -b $prefix.step4.bim -f $prefix.step4.frq -r ${hrc_sites_file} -h || true 
plink --bfile $prefix.step4 --exclude Exclude-$prefix.step4-HRC.txt --make-bed --out $prefix.step5 || true 
plink --bfile $prefix.step5 --update-chr Chromosome-$prefix.step4-HRC.txt --make-bed --out $prefix.step6 || true 
plink --bfile $prefix.step6 --update-map Position-$prefix.step4-HRC.txt --make-bed --out $prefix.step7 || true 
plink --bfile $prefix.step7 --flip Strand-Flip-$prefix.step4-HRC.txt --make-bed --out $prefix.step8 || true 
plink --bfile $prefix.step8 --a2-allele Force-Allele1-$prefix.step4-HRC.txt --make-bed --out $prefix.step9 || true 
plink --bfile $prefix.step9 --update-name ID-$prefix.step4-HRC.txt --real-ref-alleles --make-bed --out $prefix.step10 || true

**Tags:** nextflow
> **Comments:**
> **username_1**: When you say your script ran without error but you got no output files, are you looking in the workDir or publishDir? When I'm developing, I usually only write the output block declaration after I've tested the tool/script a couple of times. This is because your task process can still run successfully without an output block declaration. I can then look in the workdir to see what files have been created, etc, then finally declare them in the output block. Then I start on the next (downstream) process(es) and repeat.  
> **username_0**: The publishDir is empty and in the workDir there are only log files (chr.log, chr.step1.log, chr.step2.log etc.). It is strange!  
> **username_1**: So only files declared in your output block will get published into your publishDir. Ignore output in the publishDir while testing and developing. Or better, comment out the publishDir directive until you get your script block working in some capacity. In any case, your job script will show you what's going on. Happy to take a look if you need help.  
> **username_0**: Thank you `@Steve`. I added the job file as a png (not sure how to add text file). In line 21 `hce-ec-eur-fin-pass-qc` is the only input which create chr(1-22).(bed,bim,fam) files (for 22 chromosomes) which are used in the subsequent commands. I tried to replace `hce-ec-eur-fin-pass-qc` with `$1` but got no output. Maybe you could give a better idea.  
> **username_1**: Please correct me if I'm wrong, but it looks like this script kicks off a job array. This isn't a good fit inclusion "as-is" into a Nextflow workflow, since it won't be portable at all. I think breaking it up into two scripts is the best option here. I'll see if I can give an example of this. Also, you can paste text in a code-block. I've updated your question with the text. Click [edit](https://bioinformatics.stackexchange.com/posts/15633/edit) to see the changes I've made.  
> **username_1**: See also: https://bioinformatics.stackexchange.com/editing-help  
> **username_1**: Just realized my edit hasn't been approved yet, but you should be able to see it.  
> **username_0**: Thanks `@Steve` . Yes, it starts a job array 1-22 for each chromosome. The initial input is one single bfile (bed,bim,fam), in step1a it creates 22 bfiles (bed,bim,fam) for each chromosome. The rest of steps operate on those 22 files to create final output. Any relatable example to break the job into two would be great.  
> **username_1**: No worries at all @zillurrahman. I've updated my answer with an example on how to do this. I also refactored some of your code, hopefully it still works...  

---

## 1 answers

### Answer by username_1 on 2021-03-28. Score: 2 | Accepted answer

Using Nextflow to submit jobs this way makes little sense, since:

sbatch exits immediately after the script is successfully transferred
to the Slurm controller and assigned a Slurm job ID. The batch script
is not necessarily granted resources immediately, it may sit in the
queue of pending jobs for some time before its required resources
become available.

Assuming your job script is basically a shell script, make it executable (e.g. chmod +x PRE_PHASING_QC.JOB) and move it into a folder called bin in the root directory of your project repository. Nextflow automatically adds this folder to your $PATH, making scripts in this folder accessible in your pipeline without needing localize them and without needing to specify an absolute path to invoke them. Your code then becomes:
process prephasing {

    """
    PRE_PHASING_QC.JOB [options] 
    """
}

Finally, add process.executor = 'slurm' to your nextflow.config file. You may need to configure some other options, like resource requests, but that should get your script integrated into your pipeline without much refactoring, if any.


Contents of PRE_PHASING_QC.sh:
#!/bin/bash

set -euo pipefail

hrc_sites_file="${1}"
input_prefix="${2}"
chrom="chr${3}"

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs)

plink --bfile "${input_prefix}" --chr "${3}" --out "${chrom}" --make-bed
plink --bfile "${chrom}" --list-duplicate-vars ids-only suppress first --out "${chrom}.step1a"
plink --bfile "${chrom}" --exclude "${chrom}.step1a.dupvar" --make-bed --out "${chrom}.step1a"
plink --bfile "${chrom}.step1a" --snps-only just-acgt --make-bed --out "${chrom}.step1b"

awkcmd='
($5 == "T" && $6 == "A") || ($5 == "A" && $6 == "T") ||
($5 == "C" && $6 == "G") || ($5 == "G" && $6 == "C")
{
    print $2, "ambig"
    next
}
{
    print $2
}'

awk "${awkcmd}" "${chrom}.step1b" | grep "ambig" > "${chrom}.step1b.snplist.txt"
plink --bfile "${chrom}.step1b" --exclude "${chrom}.step1b.snplist.txt" --make-bed --out "${chrom}.step1c"

# Step 2. Zero out Mendelian errors
plink --bfile "${chrom}.step1c" --me 1 1 --set-me-missing --make-bed --out "${chrom}.step2"

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk.
plink --bfile "${chrom}.step2" --hwe 1e-6 --geno 0.02 --make-bed --out "${chrom}.step3"
plink --bfile "${chrom}.step3" --freq --make-bed --out "${chrom}.step4"

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s)
HRC-1000G-check-bim.pl -b "${chrom}.step4.bim" -f "${chrom}.step4.frq" -r "${hrc_sites_file}" -h
plink --bfile "${chrom}.step4" --exclude "Exclude-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step5"
plink --bfile "${chrom}.step5" --update-chr "Chromosome-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step6"
plink --bfile "${chrom}.step6" --update-map "Position-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step7"
plink --bfile "${chrom}.step7" --flip "Strand-Flip-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step8"
plink --bfile "${chrom}.step8" --a2-allele "Force-Allele1-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step9"
plink --bfile "${chrom}.step9" --update-name "ID-${chrom}.step4-HRC.txt" --real-ref-alleles --make-bed --out "${chrom}.step10"


Make PRE_PHASING_QC.sh executable (chmod +x PRE_PHASING_QC.sh), and replace your job script with the following. This file won't be needed by your Nextflow process, but shows you how the above shell script can be invoked:
#!/bin/bash
#SBATCH --job-name=prephase
#SBATCH --account=p697_norment
#SBATCH --time=1:00:00
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=7600M
#SBATCH --cpus-per-task=7
#SBATCH --array=1-22

set -e

. /cluster/bin/jobsetup

export PATH="/cluster/projects/p697/github/norment/moba_qc_imputation/software:${PATH}"

module load plink

cd /path/to/hce-ec-eur-fin-pass-qc

/path/to/PRE_PHASING_QC.sh \
    /cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz \
    hce-ec-eur-fin-pass-qc \
    "${SLURM_ARRAY_TASK_ID}"

Copy PRE_PHASING_QC.sh and the Oxford Perl script, HRC-1000G-check-bim.pl, into your project's bin/ directory. You may also need to make the Perl script executable. Then your Nextflow process becomes:
hrc_sites_file = file(params.hrc_sites_file)

Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size:3 )
    .set { plink_data }

Channel
    .of(1..22)
    .set { chroms }


process prephasing {

    input:
    tuple val(prefix), path(bfiles) from plink_data
    each chrom from chroms
    path hrc_sites_file

    """
    PRE_PHASING_QC.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
    """
}

Once the above runs successfully, you can declare your process outputs. Please note that none of the above was tested; hopefully not too many typos, but YMMV.

> **Comments:**
> **username_0**: Thanks `@Steve` . I think you understand my explanation and your answer makes sense. I have updated the question with new code, error and expected output. Hope it helps to debug  
> **username_1**: @zillurrahman: You get the error `Missing value declared as output parameter: chr` because you've defined 'chr' in your output block, but this value wasn't defined in your input block (or globally in your script block). You didn't show what your input block looks like. I suspect it should look something like: `tuple val(chr), val(num), val(steps), path(bfiles)`. I don't know what you job script does, but if your job script references absolute paths to files, you'll need to refactor it so that they can be passed in as variables on the command line.  
> **username_0**: My job script takes bfiles as input (e.g `plink --bfiles input_bfiles --other argument --output`). I just tried to modify this as : `plink --bfiles $1` and tried to use input channel from nextflow. `"""  PRE_PHASING_QC.job "${bfiles}" """"` In this case nothing happens, no error or success.  
> **username_1**: From our previous discussions, I thought plink took a prefix for it's '--bfile' option? Note this isn't the same as '--bfiles'. If I'm correct, I think your first task will be to work out how to supply a prefix on the command line so that it gets passed down to your plink command. In your NF process, I think you'll want: `PRE_PHASING_QC.job "${prefix}"`. And in your script I think you'll want: `plink --bfile "$1"`  
> **username_0**: Thanks `@Steve` . I tried `nextflow run prephasing1.nf --hrc_sites_file "path/to/hrc_sites.tab.gz" "path/to/hce-ec-eur-fin-pass-qc" "1"` It give me `Missing 'bind' declaration in input parameter` . Should we use `tuple(val),val(input_prefix),path(bfiles) from some_channel` ? How could I declare input channel in this particular case? It seems I could run chromosome by chromosome? That if I put "1" in nf run parameter, it will only run for chromosome 1? Should we declare it in `in_prefix=file(params.input_prefix)` and `params.chr=1` then we could change/call the arguments in command line?  
> **username_0**: Thanks `@Steve` So far I declared input_prefix and chrom as hrc_sites_file was declared in your script and specify them in command line. I got output for chromosome 1, chr1.step1a, chr1.step1b but in the next step it shows that no such file as chr1.step1b but it exists in my work/ directory. Any thought? When I use chrom 2, I got chr2.step1a, chr2.step1b and same error in the next step.  
> **username_1**: No worries @zillurrahman. I just updated the NF code to process all chromosomes like what we have done previously. I think you'd rather avoid running each chromosome one at a time. In your first comment above, you get the `Missing 'bind' declaration in input parameter` error since you're missing parameters for the `"path/to/hce-ec-eur-fin-pass-qc"` and `"1"` options.  
> **username_0**: Ohh! There was some typo in the PRE_PHASING_QC.sh script. I ran without the slurm, identfied those typos, edited them and now everything is Ok. Now, How could I declare the output files so that I get all the outputs in a single directory? So far, I tried `tuple val("${chrom}"), val(steps), path("${chrom}.${steps}.{bed,bim,fam}") into qc1_bfile` and `tuple val(prefix), val("${chrom}"), path("${prefix}-${chrom}.step4-HRC.txt") into qc1_hrc` But got error:  `Missing value decalred as output parameter: steps` How could I decalre steps so that all outputs from each step gathered in a folder?  
> **username_0**: Got it! I did as following: `tuple val("${chrom}"), path("chr${chrom}*.{bed,bim,fam}") into qc1_chr` and it works. Only thing is, file transfer taking too long. After finishing the process `Waiting files transfer to complete (2 files)` it is taking too much time. If I check output folder, I can see all my desired outputs but the waiting files transfer to complete message is still on the screen and I can not do anything in this tab. Any suggestion? However, Thanks for helping all the way!  
> **username_1**: Well done getting it working @zillurrahman. Is the issue that the file copying is taking too long, or is the issue just that while you're waiting you cannot use the current terminal? If it's the latter, it's not really a problem since the job can be run non-interactively (e.g. submitted as a job to your Slurm cluster), in the background (e.g. using screen or nohup) or can be backgrounded after i'ts started (e.g. using <ctrl-z>, then typing 'bg') to background the process.  
> **username_0**: Thanks `@Steve`. I have managed to run it using `sbatch`. It is ok now.  
> **username_0**: Hi, This process creates multiple output files `(chr1-22.step1-10.bed|bim|fam)` which are saved into qc1_chr  output channel. Now, say I need only those files which has string `step10`. How could I get that subset into a new channel? I need it into a new channel because I need to merge them with elements of other channels. Did my explanation make sense?  
> **username_1**: @zillurrahman you could use separate output channels for each step, but you might not need the other step outputs. Another way would be to select just the files you need from the qc1_chr output channel using the [filter](https://www.nextflow.io/docs/latest/operator.html#filter) operator. There's lots of ways to do this - for example: `qc1_chr.filter { bfile -> ["bed", "bim", "fam"].any { bfile.name.endsWith("step10.${it}") } }.view()`. HTH.  
> **username_0**: Thanks `@Steve`. It gives me `No such variable: Exception evaluating property 'name' for java.util.ArrayList, Reason: groovy.lang.MissingPropertyException` . If I remove `name` from the above snippet, I got `Unknown method invocation 'endsWith' on ArrayList type`. Any idea?  
> **username_1**: Sorry @zillurrahman I didn't see your qc1_chr declaration above. I think you want: `qc1_chr.filter { name, bfiles -> bfiles.collect { it.getBaseName() }.every { it.endsWith("step10") } }` HTH.  
> **username_0**: Thanks `@Steve`. It seems ok now, getting no error. Another issue. I merged other channels with it (each channel for 22 files, 1 for each chromosome) `shape_para=shape_in.merge(shape_map).merge(shape_ref).merge(shape_legend).merge(shape_sample) process phasing { input: each file(penta) from shape_para output: path("*step10*") into phased script: def (in, map, ref, legend, sample) = penta """ shapeit2 -T 32 --input-bed "${in}" --input-map "${map}" --input-ref "${ref}" "${legend}" "${sample}" --duohmm -O "${in}.phased" --seed 54321 """ }` I am getting no output or error!  
> **username_1**: I'm not sure exactly @zillurrahman, sorry. Your input declaration doesn't look right, but I'd still expect some output or log on stderr if your process is executed. If you can [ask a new question](https://bioinformatics.stackexchange.com/questions/ask), I'll see if I can add an answer. Please also show the contents of `shape_in`, `shape_ref`, `shape_legend` and `shape_sample`. Cheers.  
> **username_0**: Thanks `@Steve` I posted a new question with code for each channel here; https://bioinformatics.stackexchange.com/questions/15680/merge-multiple-channels-each-emits-multiple-files-in-nextflow  

---

