Nextflow process for samtools sort and variant calling

Question

[antoine.fauchois92]

Hi everybody !

I'm currently try to learn nextflow scripting. I'm a novice and I just begin with a script in order to sort and convert a SAM file to a BAM file and then call variants thanks to a reference sequence. So I wrote in a .nf file, the code bellow:

#! /usr/bin/env nextflow

params.ref_seq = "/home/virologie/HHV8/Ref_seq/HHV8_ref.fasta"
params.threads = 10

process sort_reads
{
    input:
    path patient

    output:
    stdout

    exec:
    println "##### Sorting mapped reads #####"

    script:
    """
    samtools view -bS ${patient} | samtools sort
    """
}

process variant_calling
{
    input:
    val bam_data

    output:
    path "/home/virologie/test.vcf.gz"

    exec:
    println "##### Calling variant #####"

    script:
    """
    bcftools mpileup -Ou -f ${params.ref_seq} ${bam_data} | \
    bcftools call -mv -Oz --threads ${params.threads} -o /home/virologie/test.vcf.gz
    """
}

workflow
{
    def patient_path = Channel.fromPath("/home/virologie/Developpement/nextflow_test/SAM/1G_S15.sam")
    bam_data = sort_reads(patient_path)
    variant_calling(bam_data)
}

The code seems to be correct, however, I get this error:

Error executing process > 'sort_reads (1)'

Caused by: Input length = 1

I don't understand why I get this and how to fix it ?

Can maybe someone help me ?

Best regards, Antoine

Answer 1

[colindaven]

You can additionally replace

samtools view -bS ${patient} | samtools sort

with something like

 samtools sort -@ params.threads -o out.sorted.bam ${patient}

Also, I'd try to not use def

def patient_path = Channel.fromPath("/home/virologie/Developpement/nextflow_test/SAM/1G_S15.sam")

patient_path = Channel.fromPath("/home/virologie/Developpement/nextflow_test/SAM/1G_S15.sam")

I would also replace the patient variable with input for clarity, and use a nextflow.config file where the input SAM and threads are defined.

[antoine.fauchois92]

Thank you ! It's working now !

Answer 2

[ATpoint]

I do not know why, but I can turn this error off by directing samtools sort to a file. I anyway recommend against using the stdout redirection you're using. This is not Unix stdout as one might assume. Everything that goes to Nextflow stdout will collected in the log files of this process, so you're basically writing large (and binary) data to these logs. If you want native Unix pipes then write the sender and receiver tools into the same process.

Essentially, use samtools sort -o someOutput.bam ${patient} -- and declare the output in output:. No need for view as sort can take SAM right away. Generally, write alignments directly to bam, like aligner (...) | samtools view -o out.bam -- that saves space and effort.