Error in using vcfutil and bcftools

[bio_informatics]

outPut.sam is my sam file from bowtie2.

I am trying to generate a consensus fq/fasta file from it.

I ran:

samtools-1.1/./samtools view -bS outPut.sam > bamOut.bam --> get the bam file out of sam file

samtools-1.1/./samtools sort bamOut.bam sorted_bam.bam --> sort the bam file

samtools-1.1/./samtools index sorted_bam.bam.bam --> create the index of the bam file sorted

Now when I try to use this command:

samtools-1.1/./samtools mpileup -uf sequence.fasta sorted_bam.bam.bam | bcftools/./bcftools view -cg - | bcftools/./vcfutils.pl vcf2fq > cns.fq

I get an error,

Error: Could not parse --min-ac g
[mpileup] 1 samples in 1 input files
Use of uninitialized value $l in numeric lt (<) at bcftools/./vcfutils.pl line 566.
Use of uninitialized value $l in numeric lt (<) at bcftools/./vcfutils.pl line 566.
<mpileup> Set max per-file depth to 8000

And everything stops.
It is something with bcftools.
References:
http://seqanswers.com/forums/showthread.php?t=11536

https://www.biostars.org/p/63429/

Please guide.

[GenoMax]

See this recent thread: http://seqanswers.com/forums/showthread.php?t=47210 If the following does not work then we can try splitting the three steps as indicated in this thread.

Can you try running your command like this:

$ samtools-1.1/samtools mpileup -uf sequence.fasta sorted_bam.bam.bam | bcftools/bcftools view -cg - | bcftools/vcfutils.pl vcf2fq > cns.fq

[bio_informatics]

Hi GenoMax,
Thanks for your response.

I ran the command you provided. It didn't work, then I proceeded with individual.

#

samtools-1.1/./samtools mpileup -uf sequence.fasta sorted_bam.bam.bam > bcf_file.bcf

ran successfully.

Below one threw error

bcftools/./bcftools view -cg bcf_file.bcf > bcftoolsout

Error is Error: Could not parse --min-ac g

When I removed -gc, it ran successfully. Strange?

bcftools/./bcftools view bcf_file.bcf > bcftoolsout

Originally Posted by GenoMax

See this recent thread: http://seqanswers.com/forums/showthread.php?t=47210 If the following does not work then we can try splitting the three steps as indicated in this thread.

Can you try running your command like this:

$ samtools-1.1/samtools mpileup -uf sequence.fasta sorted_bam.bam.bam | bcftools/bcftools view -cg - | bcftools/vcfutils.pl vcf2fq > cns.fq

[GenoMax]

It looks like the new bcftools (v.1.x) expects an integer to go with -c option. Try it with that.

-c/C, --min-ac/--max-ac <int>[:<type>] minimum/maximum count for non-reference (nref), 1st alternate (alt1) or minor (minor) alleles [nref]

[bio_informatics]

Yes, I was going through the -help.
I downloaded bcftools today.

I am constructing a fasta sequence [which is my goal].

Any random number would spoil my output.

Originally Posted by GenoMax

It looks like the new bcftools (v.1.x) expects an integer to go with -c option. Try it with that.

[bio_informatics]

Hi Genomax,
So, I tried few things with my file,

bcftools/./bcftools view --min-ac 0 -g "^miss" bcf_file.bcf

any number above 0 in --min-ac 0 isn't generating me a file worth moving ahead.
-g "^miss"

If I understand this correctly, I am excluding all the calls with missing genotypes.

Kindly throw light more on the genotype flag.

[MaximeOfOslo]

I think I have the same problem here, and I may have a clue of where it comes from, but none about how to solve it...

I'm using Illumina paired end data, that I assembled using Soap De Novo which gave me a file of the scafolds in a fasta format.

I then mapped this scaffold using bwa to a close reference, converted it to a .bam file using samtools (without forgetting to index the reference).

I'm now trying to generate the consensus file with this command

/usit/abel/u1/maxime/8_samtools/bin/samtools mpileup  -uf chrysanthemum_indicum_chloroplast.fa 1_file_sorted.bam | /usit/abel/u1/maxime/8_samtools/bin/bcftools call   -c - | /usit/abel/u1/maxime/8_samtools/bin/vcfutils.pl vcf2fq > ass_aln.fq

And i have the same error message :

[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
Use of uninitialized value in length at /usit/abel/u1/maxib/8_samtools/bin/vcfutils.pl line 565, <> line 114.
Use of uninitialized value in length at /usit/abel/u1/maxib/8_samtools/bin/vcfutils.pl line 565, <> line 114.

Note that bcftools view seems to have been replaced by bcftools call according to their page

From my understanding, I may have this problem because at some point, I lost the base quality encoding (in the de novo assembly).

Please find attached under this line the output of the bcftools call -c - command

https://dl.dropboxusercontent.com/u/...all_output.txt

[swbarnes2]

I'm still using samtools 0.1.19, can someone just post the portion of vcfutils that is causing the problem? It's a perl script, you can open it in any text editor. I feel like I've seen that error in older versions, but I can't pinpoint it without know where to look.

My first suggestion would be to reindex the reference with samtools faidx, make sure that happens without errors, as mpileup tends not to warn you if that index isn't right.

[SNPsaurus]

Can you switch to bcftools consensus?
https://samtools.github.io/bcftools/...html#consensus

Originally Posted by MaximeOfOslo

I think I have the same problem here, and I may have a clue of where it comes from, but none about how to solve it...

I'm using Illumina paired end data, that I assembled using Soap De Novo which gave me a file of the scafolds in a fasta format.

I then mapped this scaffold using bwa to a close reference, converted it to a .bam file using samtools (without forgetting to index the reference).

I'm now trying to generate the consensus file with this command

/usit/abel/u1/maxime/8_samtools/bin/samtools mpileup  -uf chrysanthemum_indicum_chloroplast.fa 1_file_sorted.bam | /usit/abel/u1/maxime/8_samtools/bin/bcftools call   -c - | /usit/abel/u1/maxime/8_samtools/bin/vcfutils.pl vcf2fq > ass_aln.fq

And i have the same error message :

[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
Use of uninitialized value in length at /usit/abel/u1/maxib/8_samtools/bin/vcfutils.pl line 565, <> line 114.
Use of uninitialized value in length at /usit/abel/u1/maxib/8_samtools/bin/vcfutils.pl line 565, <> line 114.

Note that bcftools view seems to have been replaced by bcftools call according to their page

From my understanding, I may have this problem because at some point, I lost the base quality encoding (in the de novo assembly).

Please find attached under this line the output of the bcftools call -c - command

https://dl.dropboxusercontent.com/u/...all_output.txt

[MaximeOfOslo]

Hi,
I'd like to try to use consensus. Till the faidx command, everything is fine, however, I'm now stuck at the mpileup stage.

Here is my code :

/usit/abel/u1/maxib/8_samtools/bin/samtools  view -bS 1_align.sam | /usit/abel/u1/maxib/8_samtools/bin/samtools sort - file_sorted
/usit/abel/u1/maxib/8_samtools/bin/samtools index file_sorted.bam
/usit/abel/u1/maxib/8_samtools/bin/samtools faidx chrysanthemum_indicum_chloroplast.fasta
/usit/abel/u1/maxib/8_samtools/bin/samtools mpileup -v -f chrysanthemum_indicum_chloroplast.fasta file_sorted.bam -o file.vcf.gz 
/usit/abel/u1/maxib/8_samtools/bin/bcftools consensus -i chrysanthemum_indicum_chloroplast.fasta file.vcf.gz >  out.fa

And here is the output :

[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
./sam.sh: line 4: 24661 Abandon                 /usit/abel/u1/maxib/8_samtools/bin/samtools mpileup -v -f chrysanthemum_indicum_chloroplast.fasta file_sorted.bam -o file.vcf.gz

Are the fasta lines too long ?

[GenoMax]

@Maxime: Are you using new versions of samtools and bcftools (v.1.x)? You can't mix and match old/new versions.

[MaximeOfOslo]

Thanks for your help @GenoMax, it's highly appreciated
Yes, I'm using the latest versions of both samtools and bcftools

-bash-4.1$ ls 8_samtools/install/
bcftools-1.2  htslib-1.2.1  samtools-1.2

The fasta file I'm using as a reference is the coding sequence of this chloroplast genome : http://www.ncbi.nlm.nih.gov/nuccore/JN867589.1

[GenoMax]

Was the NCBI sequence used for generating the index and doing alignments? Is that the only error you are seeing?

[MaximeOfOslo]

Indeed it is. The alignement was performed with bwa using the fasta file mentioned above.

Here is the complete output of my script :

-bash-4.1$ ./sam.sh 
[bam_sort_core] merging from 20 files...
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
./sam.sh: line 4:  4949 Abandon                 /usit/abel/u1/maxib/8_samtools/bin/samtools mpileup -v -f chrysanthemum_indicum_chloroplast.fasta file_sorted.bam -o file.vcf.gz

About:   Create consensus sequence by applying VCF variants to a reference
         fasta file.
Usage:   bcftools consensus [OPTIONS] <file.vcf>
Options:
    -f, --fasta-ref <file>     reference sequence in fasta format
    -H, --haplotype <1|2>      apply variants for the given haplotype
    -i, --iupac-codes          output variants in the form of IUPAC ambiguity codes
    -m, --mask <file>          replace regions with N
    -o, --output <file>        write output to a file [standard output]
    -c, --chain <file>         write a chain file for liftover
    -s, --sample <name>        apply variants of the given sample
Examples:
   # Get the consensus for one region. The fasta header lines are then expected
   # in the form ">chr:from-to".
   samtools faidx ref.fa 8:11870-11890 | bcftools consensus in.vcf.gz > out.fa

Up until mpileup, there are no error. After samtools consensus print its error message because no vcf file was generated at the mpileup stage.

[GenoMax]

Have you tried to manually run the mpileup command? Are all the files in the current directory? Does the sorted bam file look to be of about the same size?

[MaximeOfOslo]

The sorted and not-sorted bam files are the same size

-bash-4.1$ pwd
/usit/abel/u1/maxib/1_data/1_project/1st_assembly_strategy
-bash-4.1$ du -sh *
7,0G	1_align.sam
84K	chrysanthemum_indicum_chloroplast.fasta
3,5K	chrysanthemum_indicum_chloroplast.fasta.fai
15G	contig.fa
1,8G	file.bam
1,8G	file_sorted.bam
6,5K	file_sorted.bam.bai
0	file.vcf.gz
0	out.fa
512	sam.sh
6,3G	scafseq.fa
1,5K	test.vcf.gz
0	vcffile

Running manually mpileup produces the same error

-bash-4.1$ /usit/abel/u1/maxib/8_samtools/bin/samtools mpileup  -v -f chrysanthemum_indicum_chloroplast.fasta file_sorted.bam -o file.vcf.gz 
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
Abandon

[GenoMax]

By chance do you have extremely deep coverage (> 8000)? That is a small genome and the result is a large bam file.

[MaximeOfOslo]

Well, I'm working with WGS data to extract, for the moment, the chloroplast genome.
So I have 307 210 727 reads of mean length 151 which equals 46 696 030 504 base pairs.
The chloroplast I've mapped them to is 86444 bp.
So the coverage is around 540 188...

Well, I think you found the problem ! Thanks for your help, I'll randomly subsample my fastq files before alignment by 1000 folds !
ps : to whom might be interested, here is a script to do it :

# Written by  Aaronquinlan
# https://www.biostars.org/p/6544/
# Starting FASTQ files
export FQ1=1.fq
export FQ2=2.fq

# The names of the random subsets you wish to create
export FQ1SUBSET=1.rand.fq
export FQ2SUBSET=2.rand.fq

# How many random pairs do we want?
export N=100

# paste the two FASTQ such that the 
# header, seqs, seps, and quals occur "next" to one another
  paste $FQ1 $FQ2 | \
# "linearize" the two mates into a single record.  Add a random number to the front of each line
  awk 'BEGIN{srand()}; {OFS="\t"; \
                        getline seqs; getline sep; getline quals; \
                        print rand(),$0,seqs,sep,quals}' | \
# sort by the random number
  sort -k1,1 | \
# grab the first N records
  head -n $N | \
# Convert the stream back to 2 separate FASTQ files.
  awk '{OFS="\n"; \
        print $2,$4,$6,$8 >> ENVIRON["FQ1SUBSET"]; \
        print $3,$5,$7,$9 >> ENVIRON["FQ2SUBSET"]}'

[GenoMax]

540,000x

Reformat.sh from BBMap can also subsample directly from sam/bam.

This can save you some alignment time.

$ reformat.sh in=in.sam out=out.sam sample=some_number_here

[swbarnes2]

Picard tools can also randomly downsample a .bam file.

And the cheesy way to do it yourself would be to use awk or grep to only grab reads from a particular tile, one not on the edge would be preferable.