SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
GATK: sorting vcf file given a reference file jorge Bioinformatics 4 01-14-2015 12:16 PM
What is a BCF and VCF file? What do they stand for? prs321 Bioinformatics 2 06-25-2013 12:27 PM
Problem converting bcf to vcf cowman Bioinformatics 0 03-28-2013 03:43 AM
Clarification on BCF file format phatjoe Bioinformatics 2 01-09-2012 04:22 PM
Converting Dindel VCF file to GATK BED file MolecularToast Bioinformatics 2 09-24-2011 06:38 PM

Reply
 
Thread Tools
Old 07-02-2013, 10:51 AM   #1
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default How do I go about converting a BCF file to a VCF file?

I tried using BCFtools but I keep getting an error.



1. I had 2 SAM files that converted to BAM, sorted by chromosome, and finally indexed using Picard.




2. Using the 2 manipulated BAM files, I used the mpileup command in SAMtools. Here is the specific command:

samtools mpileup -f ref.fa in.bam in2.bam > in_in2_pileup.bcf




3. After this step I wanted to get the vcf format. I used bcftools and the command view in order to do this:

bcftools view in_in2_pileup.bcf > in_in2_pileup.vcf


The error that I received after executing this command:

incorrect number of fields (0 != 5) at 0:0






PS. I tried looking my problem up in google search and every example seemed to be irrelevant in my case.
prs321 is offline   Reply With Quote
Old 07-02-2013, 11:25 AM   #2
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Every example? How about http://samtools.sourceforge.net/mpileup.shtml. One of my first hits via Google.

However what you are probably overlooking is the proper command line option to mpileup. Either look at the above web page and/or type in 'samtools mpileup' and read the help or look at your output file. It will be good for your eyes to spot the mistake. :-)
westerman is offline   Reply With Quote
Old 07-02-2013, 11:36 AM   #3
prs321
Member
 
Location: US

Join Date: Jun 2013
Posts: 96
Default

Ahh I see i see, thanks.
prs321 is offline   Reply With Quote
Old 03-05-2014, 11:17 PM   #4
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

I have the same problem but I have not "seen" yet ...

i want to use mpileup without -u and -g options but all the command line used did not work..
Giffredo is offline   Reply With Quote
Old 03-06-2014, 05:21 AM   #5
TiborNagy
Senior Member
 
Location: Budapest

Join Date: Mar 2010
Posts: 329
Default

Can you share us the command line?
TiborNagy is offline   Reply With Quote
Old 03-06-2014, 05:31 AM   #6
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

Code:
samtools mpileup -f .fas .sorted.bam > .bcf | bcftools view > .vcf&

Just now I have used

Code:
samtools mpileup -f .fas .sorted.bam | tee teeoutput.bcf&
This "artefact" maybe works..
Giffredo is offline   Reply With Quote
Old 03-06-2014, 06:51 AM   #7
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

I just replied on biostars, but this won't work since bcftools is expecting BCF, and you're giving it mpileup, which is text.
dpryan is offline   Reply With Quote
Old 03-07-2014, 03:29 AM   #8
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

So, are you telling me that the output from
Code:
samtools mpileup -f .fas .sorted.bam
is a .txt?

.. then if I use..

Code:
samtools mpileup -f .fas .sorted.bam | tee teeoutput.txt&
...I ll reach the table that I pine for!!
Giffredo is offline   Reply With Quote
Old 03-07-2014, 04:33 AM   #9
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by Giffredo View Post
.. then if I use..

Code:
samtools mpileup -f .fas .sorted.bam | tee teeoutput.txt&
...I ll reach the table that I pine for!!
Well, it depends on what sort of table you're after. If you just want the pileup and not variant calls then:

Code:
samtools mpileup -f .fas .sorted.bam > output.txt
would seem to do what you want. There'd be no need to pipe things to tee in that case. Also, there's no reason to always end commands with "&", particularly if they'll just be printing stuff to the screen.
dpryan is offline   Reply With Quote
Old 03-07-2014, 04:58 AM   #10
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

Ok thanks!!! I will try.
About & I know.. I put it because of force of habit.. ...
Giffredo is offline   Reply With Quote
Old 03-07-2014, 05:14 AM   #11
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

chrM 136 A 6 ,,,,,, 896774
chrM 137 A 6 ,,,,,, ?=@===
chrM 138 A 6 ,,,,,, ?=<=8<
chrM 139 C 6 ,,,,,, 887801
chrM 140 C 6 ,,,,,, @>==9;
chrM 141 C 6 ,,,,,, ><==9;
chrM 142 T 6 ,,,,,, 747701

the manual say: each line represents a genomic position, consisting of chromosome name, coordinate, reference base, read bases, read qualities and alignment mapping qualities. Information on match, mismatch, indel, strand, mapping quality and start and end of a read are all encoded at the read base column.

it is not like expected... this txt is not useful at all.
Giffredo is offline   Reply With Quote
Old 03-07-2014, 06:43 AM   #12
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

That text file is exactly the output that the manual is describing, so I'm not sure what else you were expecting.
dpryan is offline   Reply With Quote
Old 03-07-2014, 07:20 AM   #13
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

I think it is not...
I expected something like this but with the right symbol in the right positions in order to have the possibility to translate and understand them!
where are the read bases? and the quality value? why so many ',' if the position is one?
and..
what mean @ it is not in the legend.. for me there is characters mismatch because of file wrong conversion...

Comes on the question: using pmileup is possible to reach other types of results different from these i reached so far?
Giffredo is offline   Reply With Quote
Old 03-07-2014, 07:27 AM   #14
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

You need to read the manual a bit more. "," is a base call, it just means "the same as the reference, on the reverse orientation". The base quality scores are in the last column, they're phred encoded. "@", for example, is 31.

Regarding mpileup vs. pileup, the output is the same. The only difference is that mpileup can deal with multiple files at once (it just tacks on an extra 3 columns per file).

Perhaps it would help if you mentioned what your actual goal is.
dpryan is offline   Reply With Quote
Old 03-10-2014, 12:10 AM   #15
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

Quote:
"," is a base call, it just means "the same as the reference, on the reverse orientation". The base quality scores are in the last column, they're phred encoded. "@", for example, is 31.
OK.. but I have still some doubts: I know "," is a base call and for this reason I don t understand why I have more than one "," for only one base in one position...

..and what is Phred code? I know the ASCII code but Phred for me is only a number derived from a logarithmic operation...

My goal is measure the editing sites position, splicing, post transcription modifications in general from my mutant sample mRNA. And I want to make a statistical work on these variations.
Giffredo is offline   Reply With Quote
Old 03-10-2014, 01:26 AM   #16
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

You'll have a call for every read covering that site, note the depth column. While Phred scores are just log10 values, they're usually encoded in ASCII by first adding 33 to the Phred score. The information provided in a pileup won't be useful for splicing, but it could help with finding RNA editing and other modifications (at least if they're observable in sequencing). Whether you want to use the pileup or just the VCF is sort of up to you.
dpryan is offline   Reply With Quote
Old 03-11-2014, 03:25 AM   #17
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

ok I understood the output that mpileup give me.

Now I think that should exist a "script" to translate the output.. I mean: if I have ....,. in the output it means I have 6 reads in total in that X position in which 5 reads are forward and 1 is reverse. I think that for a normal py script, to do this kind of translation could be easy..
Giffredo is offline   Reply With Quote
Old 03-11-2014, 05:05 AM   #18
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Yeah, that'd be pretty straight-forward to do in most languages. I should note that there is a C API for this whole process if you prefer that.
dpryan is offline   Reply With Quote
Old 03-12-2014, 04:32 AM   #19
Giffredo
Member
 
Location: Paris

Join Date: Feb 2014
Posts: 36
Default

before starting to make a script, I tried VarScan..

Reading input from Pileup_merge.pileup (VARSCAN)
chrM 1 G 4 3 G:3:1:21:1:0:3:0
chrM 2 G 4 4 G:4:1:24:1:0:4:0
chrM 3 A 4 4 A:4:1:24:1:0:4:0
chrM 4 T 4 3 T:3:1:23:1:0:3:0
chrM 5 C 2 1 C:1:1:24:1:0:1:0
chrM 6 C 4 3 C:3:1:23:1:0:3:0

mpileup file:
chrM 1 G 4 ^*,^7,^7,^*, 5529
chrM 2 G 4 ,,,, 7;8=
chrM 3 A 4 ,,,, 7:7=
chrM 4 T 4 ,,,, 465=
chrM 5 C 2 ,, 19
chrM 6 C 4 ,,,, 465=

A part that VARSCAN give me error after 300 positions (I ll think after about this error), I don t see the link between these 2 output..
Giffredo is offline   Reply With Quote
Old 03-12-2014, 06:37 AM   #20
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

You might consider starting another thread. I'm unfamiliar with varscan's output and you're more likely to find someone who is if you start a thread mentioning varscan in the subject.
dpryan 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 04:36 PM.


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