SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools mpileup generates an empty output david.tamborero Bioinformatics 6 07-21-2014 08:40 AM
Breakdancer_max output is empty inou13 Bioinformatics 4 04-25-2012 03:54 PM
View TopHat output xiaohua Bioinformatics 2 04-11-2012 05:11 PM
genotype in the output file of bcftools pengchy Bioinformatics 2 03-24-2012 10:52 AM
Bfast output and "Empty Sequence Dictionary" in .sam output aiden Bioinformatics 1 05-28-2010 06:50 PM

Reply
 
Thread Tools
Old 02-07-2012, 07:46 AM   #1
manore
Member
 
Location: paris

Join Date: Jun 2011
Posts: 19
Default bcftools view generates an empty output

Hello,


file.sam
#I convert SAM to BAM
Code:
samtools view -bS file -o file.bam
#I sort BAM file
Code:
samtools sort file.bam file.srt
I remove duplicat, then I detect variants with samtools mpileup

Code:
samtools mpileup -uf genome.fasta srt_file.bam | bcftools view -bvcg - > file_raw.bcf
I obtain no errors, bcf file contain data. (file is about gigabase). => samtools mpileup generates an empty output.

When I convert bcf into vcf, I obtain
Code:
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-sample depth to 8000
and different lines with
Code:
[bcfview] 434500000 sites processed.
[afs] 0:54828.077 1:20102.090 2:25069.833
[bcfview] 434600000 sites processed.
[afs] 0:44750.569 1:24548.359 2:30701.072
but the vcf file only contains the header and this line :
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT srt_file.bam

#convert bcf into vcf
Code:
bcftools view file_raw.bcf | vcfutils.pl varFilter -D100 > variant.flt.vcf
Someone can help me?
Thanks in advance
manore is offline   Reply With Quote
Old 02-07-2012, 09:35 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Sounds like it called no SNPs. Maybe you should generate the pileup, and eyeball some of it it to see if everything looks alright.
swbarnes2 is offline   Reply With Quote
Old 02-08-2012, 03:00 AM   #3
manore
Member
 
Location: paris

Join Date: Jun 2011
Posts: 19
Default

I decide to proceed step by step,

Code:
samtools mpileup -uf genome.fasta srt_file.bam > output_mpileup
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000

ls -lah sortie_mpileup
xxx 35G Feb x xxx output_mpileup

Code:
bcftools view -bvcg output_mpileup > output_bcf
ls -lah output_bcf
xxx 1.3G Feb x xxx output_bcf
bcftools view output_bcf | more

output_bcf file is not empty!

Code:
bcftools view output_bcf | vcfutils.pl varFilter -D100 > variant.flt.vcf
variant.flt.vcf file only contains the header and this line :
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT srt_file.bam


I don't know where is the problem!
manore is offline   Reply With Quote
Old 02-08-2012, 08:29 AM   #4
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Again, why don't you manaully inspect the pileup? Is filtering everything with a coverage higher than 100 really appropriate for your dataset?
swbarnes2 is offline   Reply With Quote
Old 02-13-2012, 12:47 AM   #5
manore
Member
 
Location: paris

Join Date: Jun 2011
Posts: 19
Default

Thanks,

The problem is that vcf output indiscates only N as an SNP!

Code:
less xxx.vcf
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT C01HGABXX_002_s_7_brut_rm.bam
chromosome:xxx 400 . N A 3.02 . DP=1;AF1=1;CI95=0.5,1;DP4=0,0,1,0;MQ=35;FQ=-30 PL 30,3,0
manore is offline   Reply With Quote
Old 02-13-2012, 05:34 AM   #6
manore
Member
 
Location: paris

Join Date: Jun 2011
Posts: 19
Default

My problem is that the pileup file contains all "N's" in the reference sequence field (column 3).
I am mapping reads to a genome reference using stampy. I have converted the bwa sam alignment file to a sorted bam file using samtools. When I try running pileup to identify variants, I get only N's in the reference base position

*I have checked the reference fasta file and don't find anything wrong.
* the headers are consistent between SAM and the original fasta.

see http://seqanswers.com/forums/showthread.php?t=5506

Someone can hepl me?
manore is offline   Reply With Quote
Old 02-13-2012, 08:29 AM   #7
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Double check the reference fasta index. You can make it with samtools faidx.
That index has to be made correctly for mpileup to work right. mpileup will run that command itself if it doesn't see the fasta index, and it will say something if that step fails, but it will keep going on anyway without it, and it will output all N's like you see.
swbarnes2 is offline   Reply With Quote
Old 02-15-2012, 12:12 AM   #8
manore
Member
 
Location: paris

Join Date: Jun 2011
Posts: 19
Default

I make a index reference sequence
Code:
samtools faidx genome.fasta
It creates genome.fai file. Hovewer, I obtain the same results!
manore is offline   Reply With Quote
Old 02-15-2012, 08:55 AM   #9
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Are 100% sure that the reference you indexed is the same as the file you aligned to? Spot check the names in the .sam, make sure that they match perfectly with the names in the fasta and the fasta.fai.

You definately can't call SNPs until you get mpileup to reconize the association bewteen the reference you aligned to, and the reference fasta you are pointing it to, and having N's in the reference allele indicates that it's not making that association yet. Are there any odd non-space chracters in your chromsome name that maybe your aligner is truncating when making the .sam?
swbarnes2 is offline   Reply With Quote
Old 02-19-2012, 11:56 PM   #10
manore
Member
 
Location: paris

Join Date: Jun 2011
Posts: 19
Default

I used v0.1.16 (r973) version of samtools; and now with v0.1.17 version, the pileup file doesn't contain "N" in the reference sequence field.
Older versions of SAMtools support the pileup command. As of SAMtools v0.1.17 (r973) and later, the pileup command is deprecated and has been replaced with mpileup to accommodate multi-sample calling.

Thanks for your help swbarnes2.
manore is offline   Reply With Quote
Old 03-12-2012, 04:04 AM   #11
wzhjlau2009
Junior Member
 
Location: china

Join Date: Jun 2011
Posts: 5
Question

how about the problem now?solved?
wzhjlau2009 is offline   Reply With Quote
Old 06-26-2012, 04:46 AM   #12
Anjali
Member
 
Location: USA

Join Date: Dec 2011
Posts: 17
Default

Hi,
I am also facing the similar problem.
I use samtools 1.1.18 and with mpileup command I get a bcf file but the vcf file is only with the header
when I try changing the varFilter read depth, the vcf file is of 0 size

Can any one of you please help me
Anjali is offline   Reply With Quote
Old 01-23-2016, 06:58 PM   #13
Jeffzheng
Junior Member
 
Location: CA

Join Date: Jan 2016
Posts: 1
Angry have this problem solved?

I also got empty vcf file
Jeffzheng is offline   Reply With Quote
Old 02-05-2016, 06:34 AM   #14
dvdhover
Junior Member
 
Location: Georgia

Join Date: Oct 2015
Posts: 1
Default htslib matters

I got into the same problem.
Code:
samtools mpileup -ugf <ref.fa> -Q 0 <bamfile.bam> | bcftools call -vc | less
I can see the called variants.
Code:
samtools mpileup -ugf <ref.fa> -Q 0 <bamfile.bam> | bcftools call -vc > <mutation.vcf>
And I get an empty vcf file. As have been pointed out by others, the problem lies in exporting the vcf file. This problem almost drove me crazy. I reinstalled the samtools and bcftools. Still not working. Finally, I installed htslib. It worked! It turned out that htslib is not optional! So you may need to make sure htslib is installed in your machine. Good luck.
dvdhover is offline   Reply With Quote
Old 07-08-2016, 02:19 AM   #15
mmmm
Senior Member
 
Location: UK

Join Date: Jul 2013
Posts: 131
Default bcftools view -b

Dear all, what is wrong in the bcftools view -bvcg ? (why -b option can not be found?). thanks
mmmm is offline   Reply With Quote
Reply

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 01:30 AM.


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