Workflow for Nanopore Sequencing of 10x single cell libraries

public 1yr ago Version: v1.1.0 0 bookmarks

View Workflow

ScNapBar (single-cell Nanopore barcode demultiplexer) is a workflow to assign barcodes to long-read single-cell sequencing data. ScNapBar enables cell barcode assignment with high accuracy using unique molecular identifiers (UMI) or a Naïve Bayes probabilistic approach. It requires bam files from both Nanopore and Illumina reads, then builds a model based on the parameters estimated from the two libraries.

If you use ScNapBar , please cite the following paper:

Wang Q, Boenigk S, Boehm V, Gehring NH, Altmueller J, Dieterich C. Single cell transcriptome sequencing on the Nanopore platform with ScNapBar. RNA. 2021 Apr 27;27(7):763–70. doi:10.1261/rna.078154.120 .

Installation

Software dependencies are managed using conda , for more information see
https://docs.conda.io/projects/conda/en/latest/user-guide/install/ .

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
conda config --set auto_activate_base false

Install scNapBar .

# Clone repository, add --recursive to include the sequan submodule.
git clone --recursive https://github.com/dieterich-lab/single-cell-nanopore.git
cd single-cell-nanopore
# Create environment...
conda env create --name scNapBar --file environment.yaml
conda activate scNapBar
cmake .
make

Quick run

Test scNapBar using the example data under data . Download the reference genome and place the file under data . Edit the provided config.yaml file by changing the paths to output and temporary directories, as well as the name of the reference genome. Run the snakemake command under the conda environment. Use the -j parameter to specify the number of available cores.

snakemake -j 12 --printshellcmds --verbose

You can also submit the job via job schedulers. We have provided an example using SLURM. Adjust the cluster.json file, or use your own snakemake profile .

snakemake -j 12 --until run_umi_seq --printshellcmds --verbose --cluster-config cluster.json --cluster "sbatch -A {cluster.account} --mem={cluster.mem} -t {cluster.time} -c {cluster.threads} -p {cluster.partition}"

scNapBar general usage. Edit the provided config.yaml file to match your own sequence files, reference genome, annotations, etc . Update the adapter and polyT length that fit your libraries. Run the snakemake command under the conda environment.

Two modes for the cell barcode assignment

ScNapBar (option 1 default) uses a probabilistic model for barcode assignment, which performs very well in cases of low sequencing saturation.
ScNapBar (option 2) assigns barcode based on matched Illumina UMIs without additional probabilistic modeling. Use the following command to start this mode:

snakemake --until run_umi_seq

Input files:

For testing purposes, we suggest downloading GSE130708 or PRJNA722142 ( we used chr17 as an example dataset in the data folder, see Quick run ).

We use the following naming convention:

possorted_genome_bam.bam from Cell Ranger, or barcodes.tsv.gz , features.tsv.gz , and matrix.mtx.gz from the filtered_feature_bc_matrix folder of Cell Ranger.
Background cell barcodes re-named as barcodes_raw.tsv.gz , from the barcodes.tsv.gz of the raw_feature_bc_matrix folder of Cell Ranger.
Nanopore reads ( Nanopore.fq.gz ) in FASTQ (compressed) format.
Reference genome in FASTA format.
Annotation file in refFlat format (for option 2). Note: The file must have the .refFlat extension. If running on the example data, the file must first be uncompressed. See also the UCSC genome annotation database .

Output files:

The output files are written to the results folder (this directory must exist before running the pipeline), or any folder specified by dir_out in the config.yaml file. Example output files from the quick run using the example data are provided under analysis . We also included SRSF2 isoforms characterized by barcodes from the published manuscript (SRSF2.GFP+.bam and SRSF2.GFP-.bam).

The main target output files are real.label (option 1) or real.umi (option 2):

real.label : read_id, barcode, score. The barcode assignment from the real Nanopore reads with scores. The scores range from the cutoff set in config.yaml to 99. Reads assigned to multiple barcodes are removed if both are above the score cutoff.
real.umi : The barcodes assignment of the Nanopore reads with matched Illumina UMIs from the same cell and the same gene.

Other files include:

Nanopore.bam : The mapping of the real Nanopore reads generated by minimap2. Note : the name of this file is given by the basename value of the config key nanopore_fq .
sim.prob and real.prob : The complete feature tables used to generate the probability scores of the simulated reads and the real Nanopore reads, respectively (option 1).
sim.label : the barcode assignment from the simulated reads with scores. The scores range from 0-99, and larger scores indicate higher confidence for the assignment. Reads assigned to multiple barcodes only have the assignment with the highest score retained (option 1).
sim_barcodes.txt : The ground-truth of cell barcodes in the simulated reads (option 1).
sim.model.rda : The naive bayes model trained from the simulated reads (option 1).
genome.fa : The artificial genome used for generating the simulated reads (option 1).
sim_umi.fasta and real_umi.fasta : The rest DNA sequences after removing the adapter and barcode sequences.

Parameters:

Parameters of `ScNapBar`

The parameters are set in the config.yaml . If the entry is a file, then it must be placed under the data folder.

reference_genome : The reference genome sequences in FASTA format.
nanopore_fq : Nanopore reads you want to process in FASTQ (compressed) format.
adapter : 10x genomics P1 adapter sequences.
polyTlength : Number of poly-Ts you want to simulate.
cdnalength : Number of nucleotides used to append to each entry in the artifical genome.
umilength : Number of nucleotides of the UMI sequences.
barcodelength : Number of nucleotides of the cell barcode sequences.
numSimReads : Number of Nanopore reads to simulate.
numSimReads : Number of Illumina reads to sample from the Illumina sequencing.
cutoff : Score cutoff of the barcode assignment of the real Nanopore reads.
percent_raw : A fraction number representing the percentage of additional simulated reads you want to use as true negatives. These reads contain the cell barcodes from the background rather than the whitelist. From our experience, there are about 20% reads do not contain the cell barcodes from the whitelist in the 10x genomics library.
threads : Number of CPUs for the multiple-threaded jobs. Note : See Control the number of cores/threads per rule
cdnaseq : DNA sequences used to append to each entry in the artificial genome.
nano_seed : Seed for the pseudo-random number generator (NanoSim).

Parameters of `singleCellPipe`

singleCellPipe is modified from flexbar which performs the adapter and barcode alignments. Most of the flexbar parameters are available in this program, yet there are a few distinctive parameters as follows:

-ul : Number of nucleotides of barcode and UMI sequences. E.g., the parameter should be 26 for a 16bp barcode and 10bp UMI library.
-kb : Number of additional nucleotides to search after the adapter alignment.
-fl : Number of nucleotides from both ends for searching the barcode sequences.

Step-by-step instructions

In this paragraph, we explain the use of each major job in the pipeline.

Parameter estimation from Illumina data

build_illumina : If Illumina bam file is not provided, we use the filtered cellular barcodes ( barcodes.tsv.gz ) detected from Cell Ranger pipeline, and produce an Illumina bam file based on the frequencies in the matrices.
find_dist : add some more background cellular barcodes into the cell barcode whitelist, to make sure we do not align the read to sub-optimal barcodes due to the absence of the real barcode sequences. It retrieves all the other barcodes within two edit-distances from the filtered cellular barcodes.
get_cbfreq : use the read counts for each barcode as prior knowledge in the Bayesian model.
align_longreads : We align the Nanopore reads to the reference genome using Minimap2 with long-reads settings .

Build predictive model

build_genome : We generated an artificial "genome" which contains only the cDNA primer from 10x Chromium Single Cell V3 , cellular barcode and UMI sequences as the same counts as the Illumina library, followed by 20bp oligo-dT and 32bp cDNA sequences in our pipeline, in order to estimate the likelihood of barcode mismatches and indels in our model.
build_nanosim : The Nanopore error profile is produced using the "read_analysis.py" from NanoSim, and creates a directory under the analysis folder that contains the error profile of the Nanopore reads.
sim_reads : We generate a number of Nanopore reads based on the artificial we built previously using NanoSim.
build_test : The generated Nanopore reads were trimmed to 100bp and wrote to a bam file, and the ground truth barcode sequences can be known by looking into the corresponding genomic locations from the artificial genome in the pipeline.
run_pipe_sim and run_pipe_real : run the feature extraction pipeline on the bam file of the simulated and the real Nanopore reads, respectively.
add_label : By comparing to our ground truth barcode sequences, we assign either 0 or 1 as labels to each barcode alignment, indicating whether the corresponding alignment is correct or not.
build_model : build a naive bayesian model based on the label and previously extracted barcode alignment features.
pred : Use the naive bayesian model previously built and predict the likelihood given the alignment features from the other Nanopore reads sequenced from the same cDNA library. Then we use the Bayesian theorem to calculate the posterior probabilities that the barcode alignment is correct among all potential barcodes. The predicted probabilities allow benchmarking our predictions with the other simulated reads, or do barcode assignment with the real Nanopore reads.
get_gene_umi : create UMI whitelist for each gene based on the Illumina data.
run_umi_sim : perform UMI alignment against the UMI whitelist of the same gene on the sequences that have the barcode trimmed.
filter_sim and filter_pred : output the reads which passed the cutoff of the predicted probabilities.

Installation troubleshooting

Q: fatal error: tbb/pipeline.h: No such file or directory when compiling singleCellPipe .

A: Please run conda install [--name scNapBar -c conda-forge] tbb=2020.3 tbb-devel=2020.3 to install the required TBB libraries.

Q: fatal error: seqan/basic.h: No such file or directory when compiling singleCellPipe .

A: Please download SeqAn first and move the SeqAn include folder to seqan , or make sure you clone the repository with the --recursive flag. If you forgot this flag, you can always git submodule update --init afterwards.

Authors

Qi Wang < [email protected] >
Sven Bönigk < [email protected] >
Christoph Dieterich

Currently maintained by Etienne Boileau < [email protected] >

Code Snippets

options(stringsAsFactors = FALSE)
Sys.setlocale("LC_NUMERIC","C")
args = commandArgs(trailingOnly=TRUE)
x=read.table(args[1],sep="\t",header=TRUE)
y=read.table(args[2],sep="\t",header=FALSE,row.names=1)
x[,ncol(x)] = 1-as.integer(y[x[,1],2]==x[,2])
colnames(x)[ncol(x)]='label'
write.table(x,file=args[3],sep="\t",quote=F,row.names=FALSE)

R From line 1 of pipelines/add_label.r

if(!require(caret)){
    install.packages("caret",repos="https://cloud.r-project.org")
    library(caret)
}
if(!require(e1071)){
    install.packages("e1071",repos="https://cloud.r-project.org")
    library(e1071)
}
Sys.setlocale("LC_NUMERIC","C")
options(stringsAsFactors = FALSE)
args = commandArgs(trailingOnly=TRUE)
j=c(5,7,8,9,10,11) # We use feature 3-9
x=read.table(args[1],sep="\t",header=TRUE)
d=preProcess(x[,j], method=c("center","scale","BoxCox"))
x[,j]=predict(d,x[,j])
model = naiveBayes(label~., data = x[,c(j,ncol(x))])
save(model,d,file=args[2])

R caret e1071 From line 1 of pipelines/build_model.r

if(!require(stringdist)){
    install.packages("stringdist",repos="https://cloud.r-project.org")
    library(stringdist)
}
options(stringsAsFactors = FALSE)
Sys.setlocale("LC_NUMERIC","C")
args = commandArgs(trailingOnly=TRUE)
x=read.table(args[1])
y=read.table(args[2])
y=y[seq(2,nrow(y),2),1]
r=unlist(lapply(y,function(d) x[stringdist(d,x[,2],method='dl')<2,2]))
r=unique(c(r,y))
write.table(r,file=args[3],sep="\t",quote=F,row.names=F,col.names=F)

R stringdist From line 1 of pipelines/find_dist.r

if(!require(caret)){
    install.packages("caret",repos="https://cloud.r-project.org")
    library(caret)
}
if(!require(e1071)){
    install.packages("e1071",repos="https://cloud.r-project.org")
    library(e1071)
}
Sys.setlocale("LC_NUMERIC","C")
options(stringsAsFactors = FALSE)
j=c(5,7,8,9,10,11) # We use feature 3-9
args = commandArgs(trailingOnly=TRUE)
load(args[1])
y=read.table(args[2],sep="\t",header=TRUE)
y2=y
y[,j]=predict(d,y[,j])
pred = predict(model, y[,j], "raw")[,1]
# Prior knowledge from Illumina sequencing
cnt = read.table(args[3])
n = cnt[match(y[,2],cnt[,2]),1]
n[is.na(n)]=0
# Pseudocount for the barcodes which are not in the whitelist
p=1
y1 = split(cbind(y2,n,pred),y2[,1])
# Do the bayesian calculation
m=do.call(rbind,lapply(y1,function(d){
d[,ncol(d)]=unlist(lapply(seq_len(nrow(d)),function(i){
prob=d[i,'n']*d[i,'pred']/(sum(d[,'n']*d[,'pred'])+p)
as.integer(round(min(d[i,'pred'],prob)*100))
}))
d
}))
write.table(m,file=args[4],sep="\t",quote=F,row.names=FALSE,col.names=TRUE)

R caret e1071 From line 1 of pipelines/pred.r

from scipy.io import mmread
import pandas as pd
import numpy as np
import gzip
import sys

mtx_file = gzip.open(sys.argv[3], 'r')
mat = (mmread(mtx_file))
ind = mat.nonzero()
row = ind[0]
col = ind[1]
v = len(ind[0])
n = int(sys.argv[4])
l = int(sys.argv[5])
s = '*'

name_file = gzip.open(sys.argv[1], 'rt')
cell_names = name_file.read()
if type(cell_names) is bytes: cell_names = cell_names.decode('utf-8')
cell_names = cell_names.splitlines()
cell_names = np.array(['CB:Z:'+name for name in cell_names])

feature_file = gzip.open(sys.argv[2], 'rt')
features = feature_file.read()
if type(features) is bytes: features = features.decode('utf-8')
features = features.splitlines()
features = np.array(['GN:Z:'+(feature.split('\t')[1] if (len(feature.split('\t')) > 1) else feature) for feature in features])

i = np.random.choice(np.arange(v), n, replace=True, p=mat.data/sum(mat.data))
c = np.array('ATGC', dtype='c')
d = np.random.rand(n,l)*4
umi = [a.tostring().decode("utf-8") for a in c[d.astype(int)]]
print('\t'.join(['@SQ','SN:15','LN:101991189']))
pd.DataFrame({'QNAME':list(range(n)), 'FLAG':[0]*n, 'RNAME':[15]*n, 'POS':[69453703]*n, 'MAPQ':[255]*n, 'CIGAR':['100M10I']*n, 'RNEXT':['*']*n, 'PNEXT':[0]*n, 'TLEN':[0]*n, 'SEQ':[s]*n, 'QUAL':[s]*n, 'CS':['CS:Z:*']*n, 'CQ':['CQ:Z:*']*n, 'GN':[a for a in features[row[i]]],'CB':[a for a in cell_names[col[i]]],'UB':['UB:Z:'+(a) for a in umi]}).to_csv(sys.stdout, header=False, index=False, sep="\t")

Python Pandas numpy scipy From line 1 of pipelines/simreads.py

run:
  get_most_abundant_gene(input.matrix, input.feature, output.file)

SnakeMake From line 58 of master/Snakefile

shell:
  """
  java -jar -Xmx64g TagReadWithGeneExon.jar I={input.bam} O=tmpge.bam ANNOTATIONS_FILE={input.annot} ALLOW_MULTI_GENE_READS=true USE_STRAND_INFO=true
  samtools view tmpge.bam |perl -F"\\t" -ane 'print "$F[0]\\t$1\\n" if /GE:Z:(\\w+)/' > {output}
  rm tmpge.bam
  """

SnakeMake SAMtools From line 67 of master/Snakefile

shell:
  """
  gzip -dc {input} | perl -ne 'print ">$1\\n$1\\n" if /^(\\w+)/' > {output}
  """

SnakeMake From line 79 of master/Snakefile

run:
  get_reads_per_barcode(input.matrix, input.barcode, output.file)

SnakeMake From line 90 of master/Snakefile

shell:
  """
  samtools view {input}|perl -ne 'print "$2\\t$3\\t$1\\n" if /GN:Z:(\\S+).*CB:Z:([ACGT]+).*UB:Z:([ACGT]+)/'|sort|uniq > {output}
  """

SnakeMake SAMtools From line 98 of master/Snakefile

shell:
  """
  Rscript pipelines/find_dist.r {input.reads_per_barcode} {input.fa_barcode} {output}.tmp
  sort {output}.tmp|uniq|perl -ne 'print ">$_$_"' > {output}
  rm {output}.tmp
  """

SnakeMake From line 109 of master/Snakefile

shell:
  """
  minimap2 -v1 -t {threads} -2 -ax splice --MD -ub {input.ref_genome} {input.fq} | samtools sort - -@{threads} -T {params.tmp} -o {output.bam}
  samtools index -@{threads} {output.bam}
  """

SnakeMake SAMtools Minimap2 From line 125 of master/Snakefile

shell:
  """
  read_analysis.py genome -i {input.fq} -ga {input.genome_alignment} -t {threads} -o {dir_out}nanosim_model/sim
  """

SnakeMake From line 138 of master/Snakefile

shell:
  """
  python pipelines/simreads.py {input.barcode} {input.feature} {input.matrix} {params.num} {params.umilength} {params.cdnaseq} | samtools view -bS > {output}
  """

SnakeMake SAMtools From line 154 of master/Snakefile

shell:
  """
  samtools view {input.bam} | perl -ne 'print ">",++$j,"\\n" if $i%25e5==25e5-1 or $j==0;if (/CB:Z:([ACGT]+).*UB:Z:([ACGT]+)/){{print "{params.adapter}$1$2","T"x{params.polyTlength},{params.cdnaseq},"\\n";$i++}}' > {output}
  gzip -dc {input.bk_barcode} | shuf -r -n {params.num} --random-source={input.bk_barcode} | perl -ne 'print ">chr",++$j,"\\n" if $i%25e5==25e5-1 or $j==0;if (/([ACGT]+)/){{$u="";$u.=[A,T,G,C]->[rand 4]for 1..{params.umilength};print "{params.adapter}$1$u","T"x{params.polyTlength},{params.cdnaseq},"\\n";$i++}}' >> {output}
  samtools faidx {output}
  """

SnakeMake SAMtools From line 172 of master/Snakefile

shell:
  """
  awk 'NR==FNR{{a[$0];next}}!($0 in a)' <(gzip -dc {input.barcode}) <(gzip -dc {input.raw_barcode}) | gzip > {output.bk_barcode}
  """

SnakeMake From line 185 of master/Snakefile

shell:
  """
  simulator.py genome -rg {input.fa_sim} -c {dir_out}nanosim_model/sim -o {dir_out}sim -n {params.num2} -t {threads} --seed {params.seed}
  sed -i 's/;/_/g' {dir_out}sim_aligned_reads.fasta
  perl -ne '$i++ if /^>/;print if $i<={params.num}' {dir_out}sim_aligned_reads.fasta > {output}
  """

SnakeMake From line 201 of master/Snakefile

shell:
  """
  perl -ne '$L={params.readlen};chomp;if (/^>/){{$id=$_}}else{{@t=split(/_/,$id);$s=$_;$d=$t[5];if ($t[4] eq "R"){{$s=reverse $s;$s=~tr/ATGCatgc/TACGtacg/;$d=$t[7]}}$s=substr($s,$d,$t[6]);print $id,"\\n",$t[6]>$L?substr($s,$L-$t[1]%$L,$L):$s,"\\n"}}' {input} > {output.fa}
  """

SnakeMake From line 216 of master/Snakefile

shell:
  """
  perl pipelines/fa2sam.pl {input.sim} {input.barcode} {params.readlen} {params.cdnalength} | samtools view -bS > {output}
  """

SnakeMake SAMtools From line 230 of master/Snakefile

shell:
  """
  perl -ne '$L={params.readlen};next unless /^>/;$_=substr($_,1);@t=split(/_/);$d=$t[1]+$L-$t[1]%$L;print "$t[0]\\t$d\\t",$d+$L,"\\t$_"' {input.sim}|sort -k1,1 -k2,2n > {input.sim}.bed
  bedtools getfasta -fi {input.fa_sim} -bed {input.sim}.bed -name > {input.sim}.fa
  samtools view {input.bam} | perl pipelines/gtruth.pl {input.sim}.fa {params.adapterlength} {params.barcodelength} {params.umilength} {input.gene} > {output.barcode}
  cut -f1-2 {output.barcode} > {output.genes}
  """

SnakeMake SAMtools BEDTools From line 249 of master/Snakefile

shell:
  """
  if [ -f {params.prefix}_umi.fasta ]; then rm {params.prefix}_umi.fasta; fi
  bin/singleCellPipe -n {threads} -r {input.bam} -t {params.prefix} -w {input.barcode} -as {params.adapter} -ao 10 -ae 0.3 -ag -2 -hr T -hi 10 -he 0.3 -bo 5 -be 0.2 -bg -2 -ul {params.barumilength} -kb {params.keepbp} -fl 100
  awk '$2!="NA" || NR==1' {params.prefix}.tab > {output.tab}
  perl -ne 'if(/^>/){{print}}else{{chomp;print substr($_,{params.keepbp},{params.len}),"\\n"}}' {params.prefix}_umi.fasta > {output.umi}
  rm {params.prefix}.tab {params.prefix}.fasta {params.prefix}_umi.fasta {params.prefix}_paramLog.log
  """

SnakeMake From line 271 of master/Snakefile

shell:
  """
  if [ -f {params.prefix}_umi.fasta ]; then rm {params.prefix}_umi.fasta; fi
  bin/singleCellPipe -n {threads} -r {input.bam} -t {params.prefix} -w {input.barcode} -as {params.adapter} -ao 10 -ae 0.3 -ag -2 -hr T -hi 10 -he 0.3 -bo 5 -be 0.2 -bg -2 -ul {params.barumilength} -kb {params.keepbp} -fl 100
  awk '$2!="NA" || NR==1' {params.prefix}.tab > {output.tab}
  printf "Aligned to genome\t" > {output.log}
  cut -f1 {params.prefix}.tab|perl -npe 's/_end[1|2]//'|sort|uniq|wc -l >> {output.log}
  printf "Aligned to adapter\t" >> {output.log}
  awk '$3=="no"' {params.prefix}.tab|cut -f1|perl -npe 's/_end[1|2]//'|sort|uniq|wc -l >> {output.log}
  printf "Aligned to barcode\t" >> {output.log}
  awk '$2!="NA"' {params.prefix}.tab|cut -f1|perl -npe 's/_end[1|2]//'|sort|uniq|wc -l >> {output.log}
  perl -ne 'if(/^>/){{print}}else{{chomp;print substr($_,{params.keepbp},{params.len}),"\\n"}}' {params.prefix}_umi.fasta > {output.umi}
  rm {params.prefix}.tab {params.prefix}.fasta {params.prefix}_umi.fasta {params.prefix}_paramLog.log
  """

SnakeMake From line 295 of master/Snakefile

shell:
  """
  perl pipelines/get_gene_umi.pl {input.umi} {input.ge} > {output}
  """

SnakeMake From line 316 of master/Snakefile

shell:
  """
  mkdir {params.tmpdir}
  perl pipelines/umialign.pl {input.illu} {input.nano} x {params.umilength} {params.tmpdir} > {output}
  rm -fr {params.tmpdir}
  """

SnakeMake From line 330 of master/Snakefile

shell:
  """
  Rscript pipelines/add_label.r {input.tab} {input.barcode} {output}
  """

SnakeMake From line 343 of master/Snakefile

shell:
  """
  Rscript pipelines/build_model.r {input} {output}
  """

SnakeMake From line 353 of master/Snakefile

shell:
  """
  Rscript pipelines/pred.r {input.model} {input.tab} {input.reads_per_barcode} {output}
  """

SnakeMake From line 365 of master/Snakefile

shell:
  """
  Rscript pipelines/pred.r {input.model} {input.tab} {input.reads_per_barcode} {output}
  """

SnakeMake From line 377 of master/Snakefile

shell:
  """
  awk 'NR>1' {input} | sed 's/_end[1|2]//' | cut -f1-2,15 | sort -k3,3rn | awk '!_[$1]++' > {output}
  """

SnakeMake From line 389 of master/Snakefile

shell:
  """
  awk 'NR>1 && $15>{params.cutoff}' {input} | sed 's/_end[1|2]//' | awk '{{a[$1]++;b[$1]=$0}}END{{for(i in a){{if(a[i]==1)print b[i]}}}}' | cut -f1-2,15 > {output}
  """

SnakeMake From line 401 of master/Snakefile

shell:
  """
  Rscript -e 'rmarkdown::render("report.rmd",output_file="{output}")' {input.barcode} {input.slabel} {input.sprob} {input.rlabel} {input.rprob} {input.rlog} {params.cutoff} null
  """