SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
error when running samtools mpileup and bcftools linp5 Bioinformatics 5 08-22-2014 05:41 AM
samtools bcftools no FQ value StefanP Bioinformatics 2 06-20-2014 03:51 PM
BCFTOOLS 0.2 Filter problems simono101 Bioinformatics 0 05-27-2014 06:13 AM
question about bcftools merge giror Genomic Resequencing 2 02-27-2014 02:38 AM
bcftools error Morpse Bioinformatics 8 10-05-2012 12:24 PM

Reply
 
Thread Tools
Old 10-15-2014, 07:57 AM   #1
bio_informatics
Senior Member
 
Location: USA

Join Date: Nov 2013
Posts: 182
Default Error in using vcfutil and bcftools

outPut.sam is my sam file from bowtie2.

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

I ran:
Quote:
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,

Quote:

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.

Last edited by bio_informatics; 10-15-2014 at 07:59 AM.
bio_informatics is offline   Reply With Quote
Old 10-15-2014, 08:05 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

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:
Code:
$ samtools-1.1/samtools mpileup -uf sequence.fasta sorted_bam.bam.bam | bcftools/bcftools view -cg - | bcftools/vcfutils.pl vcf2fq > cns.fq
GenoMax is offline   Reply With Quote
Old 10-15-2014, 08:27 AM   #3
bio_informatics
Senior Member
 
Location: USA

Join Date: Nov 2013
Posts: 182
Default

Hi GenoMax,
Thanks for your response.

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

#
Code:
samtools-1.1/./samtools mpileup -uf sequence.fasta sorted_bam.bam.bam > bcf_file.bcf
ran successfully.

Below one threw error
Code:
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?
Quote:
bcftools/./bcftools view bcf_file.bcf > bcftoolsout
Quote:
Originally Posted by GenoMax View Post
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:
Code:
$ samtools-1.1/samtools mpileup -uf sequence.fasta sorted_bam.bam.bam | bcftools/bcftools view -cg - | bcftools/vcfutils.pl vcf2fq > cns.fq
bio_informatics is offline   Reply With Quote
Old 10-15-2014, 08:36 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

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

Quote:
-c/C, --min-ac/--max-ac <int>[:<type>] minimum/maximum count for non-reference (nref), 1st alternate (alt1) or minor (minor) alleles [nref]
GenoMax is offline   Reply With Quote
Old 10-15-2014, 08:41 AM   #5
bio_informatics
Senior Member
 
Location: USA

Join Date: Nov 2013
Posts: 182
Default

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.
Quote:
Originally Posted by GenoMax View Post
It looks like the new bcftools (v.1.x) expects an integer to go with -c option. Try it with that.
bio_informatics is offline   Reply With Quote
Old 10-15-2014, 08:52 AM   #6
bio_informatics
Senior Member
 
Location: USA

Join Date: Nov 2013
Posts: 182
Default

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

Code:
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.
bio_informatics is offline   Reply With Quote
Old 06-16-2015, 07:48 AM   #7
MaximeOfOslo
Junior Member
 
Location: Oslo

Join Date: Jun 2015
Posts: 6
Default

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

Code:
/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 :

Code:
[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 is offline   Reply With Quote
Old 06-16-2015, 11:14 AM   #8
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

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.
swbarnes2 is offline   Reply With Quote
Old 06-16-2015, 11:28 AM   #9
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 521
Default

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

Quote:
Originally Posted by MaximeOfOslo View Post
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

Code:
/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 :

Code:
[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
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Old 06-18-2015, 05:58 AM   #10
MaximeOfOslo
Junior Member
 
Location: Oslo

Join Date: Jun 2015
Posts: 6
Default

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 :
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 :

Code:
[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 ?

Last edited by MaximeOfOslo; 06-18-2015 at 06:01 AM.
MaximeOfOslo is offline   Reply With Quote
Old 06-18-2015, 06:10 AM   #11
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

@Maxime: Are you using new versions of samtools and bcftools (v.1.x)? You can't mix and match old/new versions.
GenoMax is offline   Reply With Quote
Old 06-18-2015, 06:14 AM   #12
MaximeOfOslo
Junior Member
 
Location: Oslo

Join Date: Jun 2015
Posts: 6
Default

Thanks for your help @GenoMax, it's highly appreciated
Yes, I'm using the latest versions of both samtools and bcftools
Code:
-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

Last edited by MaximeOfOslo; 06-18-2015 at 06:17 AM.
MaximeOfOslo is offline   Reply With Quote
Old 06-18-2015, 06:21 AM   #13
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

Was the NCBI sequence used for generating the index and doing alignments? Is that the only error you are seeing?
GenoMax is offline   Reply With Quote
Old 06-18-2015, 06:42 AM   #14
MaximeOfOslo
Junior Member
 
Location: Oslo

Join Date: Jun 2015
Posts: 6
Default

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

Here is the complete output of my script :
Code:
-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.
MaximeOfOslo is offline   Reply With Quote
Old 06-18-2015, 07:19 AM   #15
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

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?
GenoMax is offline   Reply With Quote
Old 06-18-2015, 07:34 AM   #16
MaximeOfOslo
Junior Member
 
Location: Oslo

Join Date: Jun 2015
Posts: 6
Default

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

Code:
-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

Code:
-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
MaximeOfOslo is offline   Reply With Quote
Old 06-18-2015, 07:57 AM   #17
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

By chance do you have extremely deep coverage (> 8000)? That is a small genome and the result is a large bam file.
GenoMax is offline   Reply With Quote
Old 06-18-2015, 08:09 AM   #18
MaximeOfOslo
Junior Member
 
Location: Oslo

Join Date: Jun 2015
Posts: 6
Default

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 :

Code:
# 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"]}'

Last edited by MaximeOfOslo; 06-18-2015 at 08:15 AM.
MaximeOfOslo is offline   Reply With Quote
Old 06-18-2015, 08:29 AM   #19
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,079
Default

540,000x

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

This can save you some alignment time.
Code:
$ reformat.sh in=in.sam out=out.sam sample=some_number_here
GenoMax is offline   Reply With Quote
Old 06-18-2015, 11:02 AM   #20
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

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.
swbarnes2 is offline   Reply With Quote
Reply

Tags
bcftools, samtools

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 11:45 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO