SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Chromosome start and end coordinates Ooinp Bioinformatics 1 04-23-2012 11:13 AM
How samtools mpileup uses paired end data to determine SNP phasing bioinformer Bioinformatics 0 03-06-2012 01:19 PM
how can i decide the exact start and end of rRNA on mitochondria bbsinfo Illumina/Solexa 0 12-02-2011 01:54 PM
Parsing a Chromosome sequence with start and end coordinates empyrean Bioinformatics 13 09-10-2011 10:23 PM
How many amounts of mRNA do you start the pair end sequencing? Chien-Yuan Chen RNA Sequencing 0 03-02-2009 10:59 AM

Reply
 
Thread Tools
Old 01-18-2013, 02:10 AM   #1
cboustred
Member
 
Location: Bristol, UK

Join Date: May 2010
Posts: 11
Default Samtools mpileup - start and end segments

Hi,

I'm having a problem calling a variant with NGS which has been confirmed via Sanger sequencing.

I'm using samtools mpileup followed by VarScan.

Bit of background
The lab are doing resequencing of human candidate genes to identify rare mutations.

Target enrichment is performed via HaloPlex, sequencing is done on a MiSeq with 150bp paired end reads.

I am aligning reads with BWA using default PE settings againts the entire human genome
(I have tried the UCSC, 1KG, and MiSeq ref genomes - all the same result)

The problem

In the pileup (generated with or without BAQ computation) there is a ~20bp region 'missing'

Example commands
samtools mpileup -f ref.fa test.bam > test_BAQ.pileup
samtools mpileup -Bf ref.fa test.bam > test_NO_BAQ.pileup

From positions chr1:76740134 to chr1:76740154 are not present in the pileup?

The bases flanking the 'missing' region are marked by ^ (start of read segment) and $ (end of read segment).

However, when I look in IGV there are plenty reads covering the 'missing' region and the variant I know is there is there clear as day! It just doesn't make it into the pileup?!

Any help / advice about what is going on here would be much appreciated

BW

Chris

Example pileup
chr1 76740129 T 51 ................................................... FF:FGBFG?FGFFEGBF6F??GBGE@FB5DF?FGDFBBFGGFD.FF?G?FG
chr1 76740130 G 51 ................................................... GF@FGFDFBFFFFEG?DBFBBGDF2DGFBFGDFGFFFDFGFFFDFF?FBFG
chr1 76740131 G 51 ...$................................................ GG9,GFDFFFGFBEGDB?FBDGFG;EGD,F>BGGBFGFFFGGDFFFBGDFG
chr1 76740132 A 50 .................................................. GG?FFDGBGFFB@G>DFFBFGFG@DGFBG?;GG?GFFFFFGFFFG>F;GF
chr1 76740133 A 50 .................................................. GGFFFFFDBFDD@G>F?FBBFF>@DGF>G??FFBFGBFGBFGFFG?>,FB
chr1 76740134 T 50 .$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$.$ >>>>>>>>>>>>9>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
chr1 76740154 A 42 ^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^],^], B=ECAACCAEBA@BEB=;C?CCBCC;>ECAECC==>>
chr1 76740155 T 42 ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, EEEGEG=GBGBD@EGEBEE;EDEGGEDEEBGDEBEGBEDDGD
chr1 76740156 A 51 ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,^],^],^],^],^],^],^],^],^], FDDGFGDGEGEG@FGEBEGDGEEGGEGEEFGEGDFGGEDFEF@BBB@BB>9
chr1 76740157 A 51 ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, FBDGFG4GEGDEEFGEDDGEGFFGGEGEEFGEGEFGGDDEGF@;DEDE@DD
chr1 76740158 T 51 ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, FBBGFG;GGG,EEFEEDBGDGFFGGEGBFEEEGEFGGF=DGF**DD9E@BD
chr1 76740159 A 51 ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, FEEGFGBGGGDDEEGE=BG;GEEDEEGFFEDEEEEEGFEEGFED@99BE>9
cboustred is offline   Reply With Quote
Old 01-21-2013, 03:56 AM   #2
cboustred
Member
 
Location: Bristol, UK

Join Date: May 2010
Posts: 11
Default

OK I have found a work around

If I use the mpileup -A flag (count anomolous read pairs) then the region is included in the pileup

BW

Chris
cboustred is offline   Reply With Quote
Old 01-30-2013, 09:27 AM   #3
dkoboldt
Member
 
Location: St. Louis

Join Date: Mar 2009
Posts: 62
Default

Chris,

I'm glad you came across the solution - it looks like there's a universal issue mapping reads in "proper" pairs at that location. Thanks for using VarScan!
dkoboldt is offline   Reply With Quote
Old 03-06-2013, 07:53 PM   #4
MWN
Junior Member
 
Location: CA

Join Date: Aug 2011
Posts: 8
Default

Quote:
Originally Posted by cboustred View Post
OK I have found a work around

If I use the mpileup -A flag (count anomolous read pairs) then the region is included in the pileup

BW

Chris
How well does Haloplex perform based on your experience? I am thinking about trying it out. However, I am worrying about FFPE DNA input requirement.
Any potential to get copy number information using Haloplex because of the multiple amplicon design? I understand the reads are not good as random reads from hybridization capture.
MWN is offline   Reply With Quote
Old 06-10-2014, 11:05 AM   #5
rpauly
Member
 
Location: Atlanta

Join Date: Apr 2011
Posts: 32
Default

I am having a similar problem with samtools mpileup:
#If I use :
samtools mpileup -r chr12:25398277-25398285 Sample1.bam >output1
where output1 is below:
chr12 25398279 N 16 G$A$*A$C$G$A$G$CCCCCCCC BB11BACBHHHHGHF0
chr12 25398280 N 9 *GGGGGGGA 1EGG?AC01
chr12 25398281 N 9 CCCCCCCCC 1EGGA/EE

These are obviously truncated lines! Why does mpileup truncate the lines in the last two columns?

Any help shall be appreciated.
~Rini
rpauly is offline   Reply With Quote
Reply

Tags
bwa, pileup, samtools mpileup, variant calling, varscan

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 03:08 PM.


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