Question: mpileup for creating a consensus sequence
0
gravatar for phernet123
19 months ago by
phernet1230
phernet1230 wrote:

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

ADD COMMENTlink modified 19 months ago by genomax67k • written 19 months ago by phernet1230
1

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"

ADD REPLYlink modified 19 months ago • written 19 months ago by genomax67k

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.

ADD REPLYlink written 19 months ago by harold.smith.tarheel4.3k

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.

ADD REPLYlink modified 19 months ago • written 19 months ago by genomax67k

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.

ADD REPLYlink modified 19 months ago • written 19 months ago by Kevin Blighe42k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 700 users visited in the last hour