SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools mpileup calls way too less SNPs TuA Bioinformatics 17 03-01-2018 05:17 PM
Impact on quality of SNP calls using samtools mpileup myi Bioinformatics 3 03-03-2014 08:18 AM
Consensus from mpileup for haploid sequences (forcing base calls - no ambiguities) ericarcher Bioinformatics 6 01-24-2014 09:52 AM
Generate consensus sequence with Samtools mpileup? Heisman Bioinformatics 4 10-28-2013 05:57 PM
Samtools mpileup calls drastically more SNPs with -I agel Bioinformatics 0 01-20-2012 02:20 PM

Reply
 
Thread Tools
Old 12-01-2012, 05:10 AM   #1
StefanP
Member
 
Location: New Zealand

Join Date: Oct 2011
Posts: 10
Default samtools mpileup consensus calls A instead of a gap

Hello,

I'm using samtools with the following command to call the consensus of my mapping:

samtools mpileup -uf Ref SEQ | bcftools view -cg - | vcfutils.pl vcf2fq > SEQ.fq

The problem is that I have sequences that have 1 and some have 2 copies of a short repeat - the reference has 2. If I use the reference with 2 copies I end up with 2 copies in the consensus since it sometimes maps the reads to the first and sometimes to the second copy of the repeat. So I exchanged one copy in the reference with N's., but end up with A's instead of the gap in the consensus.

Example:

Ref AANNNTATTATTA
Read1 AA---TATTATTA
Read2 AA---TATTATTA

it calls AAAAATATTATTA instead of AATATTATTA.

I also get the following error if I use the reference with N's instead of the second copy (that I don't get if I use the reference with 2 copies):

Unknown command "vcf2fq".


I've had the same problem with the previous pileup function!

Any idea how to change that?

Cheers,
Stefan
StefanP is offline   Reply With Quote
Old 12-01-2012, 05:41 AM   #2
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

I wonder if this is an issue with treating N in the reference as A, or down to the fact that the reads have an A before the gap? Do you have some other examples with different bases used instead of the N's?
maubp is offline   Reply With Quote
Old 12-01-2012, 06:11 AM   #3
StefanP
Member
 
Location: New Zealand

Join Date: Oct 2011
Posts: 10
Default

Hi,

Replacing "N" by "-" seams to help. It looks like it's really reading N's as A's.

Thanks!

Stefan
StefanP is offline   Reply With Quote
Old 12-01-2012, 06:16 AM   #4
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

So if you replaced the "N" with "C" or "T" or "G" do you get those letters instead of "A"? That would pretty much confirm the nature of this bug / design limitation, and give a good clue about where in the code to start looking.

Reminds me of this tweet from Vince Buffalo,
‏@vsbuffalo: Any time you only handle A, T, C, or G, the members of IUPAC kill an innocent kitten. We've all done it. #bioinfo
https://twitter.com/vsbuffalo/status/258852545634648064

Last edited by maubp; 12-01-2012 at 06:16 AM. Reason: formatting
maubp is offline   Reply With Quote
Old 12-01-2012, 06:29 AM   #5
StefanP
Member
 
Location: New Zealand

Join Date: Oct 2011
Posts: 10
Default

Hi,

Yes, I checked and it calls the reference in case the position is not covered by any read.

I should probably report that as a bug.

Thanks!

Stefan
StefanP is offline   Reply With Quote
Old 12-01-2012, 06:33 AM   #6
StefanP
Member
 
Location: New Zealand

Join Date: Oct 2011
Posts: 10
Default

It seams they changed it a little bit since the pileup function, which also interprets "-" as "A", but didn't fully resolve the problem.

Cheers,
Stefan
StefanP is offline   Reply With Quote
Old 12-01-2012, 06:55 AM   #7
StefanP
Member
 
Location: New Zealand

Join Date: Oct 2011
Posts: 10
Default

It seams to be a problem only when positions are "covered" by gaps in the read. Since it's calling N in the consensus if the position is not covered by a single read.

Stefan
StefanP 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 10:14 AM.


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