mpileup for creating a consensus sequence

Question

[phernet123]

Dear colleagues

I desperately need your help with the samtools mpileup command. I am not variant calling, rather I want to create a consensus sequence of all the mapped reads in the Sorted.bam file. I am trying to run the following bash script:

#!/bin/bash
#PBS -l select=1:ncpus=24:mem=64GB
#PBS -P CBBI0911
#PBS -q smp
#PBS -l walltime=96:00:00
#PBS -o /mnt/lustre/users/felegbeleye/my_data/stdout.femi.txt
#PBS -e /mnt/lustre/users/felegbeleye/my_data/stderr.femi.txt
#PBS -N samtools
#PBS -M --removed--
module load chpc/BIOMODULES
module load samtools/0.1.19


EXE="samtools"

ARGS="mpileup -u -f 19Dec2016_S9zPn.fasta Sorted.bam |
bcftools call -c - |
vcfutils.pl vcf2fq -D 100 > femi.psmc.fq"

cd /mnt/lustre/users/felegbeleye/my_data/Ref_Genome

${EXE} ${ARGS}

---

I am getting the following error:

mpileup: invalid option -- 'c'
open: No such file or directory
[mpileup] failed to open |: No such file or directory

I do not have an option c in my mpileup command. My indexing of the reference was done in bwa and not in Samtools. Could that be it? I feel I am doing something basic wrong but I am new at this and would appreciate some guidance.

Thank you very much for your time and help.

Femi

[GenoMax]

Try putting the entire command for ARGS on one line and see if that helps. (note: I redacted your email address from the script above).

ARGS="mpileup -u -f 19Dec2016_S9zPn.fasta Sorted.bam | bcftools call -c - | vcfutils.pl vcf2fq -D 100 > femi.psmc.fq"

[harold.smith.tarheel]

For starters, you're using an old/deprecated version (v0.1.19) that did not include the BCFtools 'call' command (I believe it was 'view' back then). Update to the current version and try again.

[GenoMax]

Good catch. phernet123 : This may just be a matter of using right module load call. Check with module avail to see if newer versions of samtools (v.1.4+) are available.

[Kevin Blighe]

My indexing of the reference was done in bwa and not in Samtools. Could that be it?

SAMtools and BWA will use different indices. If you haven't create the SAMtools indices, create it in addition to the ones for BWA with samtools faidx. Other programs use their own indices too.

I don't know if this is the exact problem. You should follow-up on what the other have already said.