1
$\begingroup$

I need to call variants on a large number of reads using both VarScan2 and BCFTools. Both of these variant calling tools use mpileup file as input to call the variants. VarScan2 command is straightforward, as we can specify the mpileup file as a command parameter. The problem is, BCFTools expects mpileup to be piped from another command and I cannot find an input parameter to specify a mpileup file that already exists (the one used with VarScan2). It seems horribly inefficient to run the same command twice instead of reusing the file obtained from previous step (especially because generating mpileup takes a substantial amount of time on my computer).

Is there a way to reuse an existing mpileup file with BCFTools?

I am using Nextflow to run my pre-processing pipeline, so could it be possible that I can pipe an input channel with mpileup into the BCFTools's call command?

I also tried doing cat file.mpileup | bcftools call ... but this does not seem to work.

Edit: cat file.mpileup | bcftools call ... does work, which would be a solution to this question, however, answers offer other, even more elegant, solutions.

$\endgroup$
2
  • $\begingroup$ I am voting to close this since the OP clarified that cat file.mpileup | bcftools does work, and there was a typo in the first attempt. $\endgroup$
    – terdon
    Sep 16 at 13:58
  • $\begingroup$ No close: There are misconceptions and nuances addressed in the answers that are of educational interest to other readers. $\endgroup$ Sep 17 at 13:17
2
$\begingroup$

Why not use tee? Please see this post.

You can write mpileup to file while also piping it to another process. For example, samtools mpileup ... | tee file.mpileup | bcftools call ...

the bcftools will run while generating the file.mpileup. After it's done, you can use file.mpileup for your VarScan2.

You can test it with the following:

touch helloworld.txt
ls | tee file.txt | grep helloworld    # output helloworld.txt
cat file.txt | grep helloworld         # output helloworld.txt
$\endgroup$
2
  • $\begingroup$ Thanks @Law for the help. It turns out cat file.mpileup | bcftools call ... worked after all. It was a silly typo mistake from my side. But still, I cannot understand how to use tee here. In your example, you have samtools mpileup ... | tee ... | bcftools call .... In my case the first part of that pipe needs to be a file instead of command. $\endgroup$
    – PPrasai
    Sep 14 at 10:28
  • $\begingroup$ In your original post, you said The problem is, BCFTools expects mpileup to be piped from another command and I cannot find an input parameter to specify a mpileup file that already exists (the one used with VarScan2). From what I understood is that bcftools expect only a piped input. You also said you didn't want to run the command twice of generating the mpileup, hence my answer. You can use tee to pipe mpileup to bcftools while also writing it to a file (aka to stdout and to file). Glad that it was only a typo then. $\endgroup$
    – user324810
    Sep 14 at 10:49
1
$\begingroup$

Bcftools and samtools both follow the very common Unix filter pattern: they read from all input files specified (as non-option plain arguments), or if there are none they read from standard input.

bcftools call is no exception:

$ bcftools call

About:   SNP/indel variant calling from VCF/BCF.
         To be used in conjunction with bcftools mpileup.
         [etc]
Usage:   bcftools call [options] <in.vcf.gz>

So you can write either of the following equivalently:

cat file.mpileup | bcftools call [options]
bcftools call [options] file.mpileup

This assumes that the file is VCF/BCF as produced by bcftools mpileup, and not actually textual “mpileup output” as suggested by the filename. If it is the latter that VarScan2 is expecting, you won't be able to share the input files between the two after all.

See this workflow example, which demonstrates bcftools call both with a pipe and a previously made input file.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.