Error while trying to get consensus Fastq from BAM
2
1
Entering edit mode
5.0 years ago

Hi all,

I am trying to get consensus Fastq from BAM file and I am running into this error. I looked for an answer in many posts but could not resolve it. Please help.

command:

/home/samtools-1.2/samtools mpileup -uf Rattus_norvegicus.Rnor_5.0.dna_sm.fa dedupRratUK_recal.bam -d 8000 | /home/bcftools-1.2/bcftools view -cg - | /home/bcftools-1.2/vcfutils.pl vcf2fq -d 10 -D 8000 > RratUK.fastq

Error:

Error: Could not parse --min-ac g
[mpileup] 1 samples in 1 input files
Use of uninitialized value in length at /home/sb47/Programs/bcftools-1.2/vcfutils.pl line 565.
Use of uninitialized value in length at /home/sb47/Programs/bcftools-1.2/vcfutils.pl line 565.

Thank you.

samtools bcftools vcfutils • 5.6k views
ADD COMMENT
0
Entering edit mode

try running step by step (make sure bam is sorted. Or try vcfutils that is attached to samtools /usr/share/samtools/vcfutils.pl)

ADD REPLY
4
Entering edit mode
5.0 years ago

the commands you're using (pileup generation + variant calling + fastq generation) seem correct, although the version of bcftools you're using doesn't work with the options you've probably seen on many old posts. the problem is that newer bcftools uses call to call variants, rather than the former bcftools view -cg

/home/samtools-1.2/samtools mpileup -uf Rattus_norvegicus.Rnor_5.0.dna_sm.fa dedupRratUK_recal.bam -d 8000 \
| /home/bcftools-1.2/bcftools call -mv \
| /home/bcftools-1.2/vcfutils.pl vcf2fq -d 10 -D 8000 \
> RratUK.fastq
ADD COMMENT
0
Entering edit mode

Just a note, the vcfutils.pl script threw an error for me. However, I was able to use this command to generate fasta consensus sequences:

samtools mpileup -uf ref.fa aln.bam | bcftools call -mv -Oz -o calls.vcf.gz
tabix calls.vcf.gz
cat ref.fa | bcftools consensus calls.vcf.gz > cns.fa

Source: https://github.com/samtools/bcftools/wiki/HOWTOs#consensus-calling

I am using samtools v1.3.1, htslib v1.3.1, and bcftools v1.3.1

ADD REPLY
0
Entering edit mode

Hi there, when I try running these commands, the line cat e1a.fasta | bcftools consensus calls.vcf.gz > EG28-SEQ-2020129-79.cns.fa gives me a Segmentation fault (core dumped). Any idea how that can be addressed?

ADD REPLY
0
Entering edit mode
5.0 years ago

fastq generation from a bam file is much easier: the reads are just there, and you just have to extract them. several tools may help you with this, but my personal choice is bedtools bamtofastq

bedtools bamtofastq [OPTIONS] -i <BAM> -fq <FASTQ>

although the newer samtools bam2fq option seems to perform fine too.

edit: this answer only applies if all the reads from the bam file are to be extracted in fastq format, but not to generate the consensus fastq requested.

ADD COMMENT
0
Entering edit mode

This is probably not what the OP is asking, the title is misleading. They probably want to get the consensus FASTA file out of the BAM file pileup.

ADD REPLY
0
Entering edit mode

you're probably right Istvan. I think it was me that I read it too fast ;)

ADD REPLY

Login before adding your answer.

Traffic: 2671 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6