SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
No output from samtool's mpileup mikheyev Bioinformatics 6 03-18-2014 05:58 AM
How to make Samtools to output both SNP and INDEL variants? shuang Bioinformatics 3 05-21-2013 11:02 AM
Using UnifiedGenotyper to call variants from haploid data? oiiio Bioinformatics 0 01-24-2012 08:29 PM
How to understand the output of mpileup like this skblazer Bioinformatics 0 12-05-2010 10:43 AM

Reply
 
Thread Tools
Old 07-19-2011, 04:37 AM   #1
gavin.oliver
Senior Member
 
Location: uk

Join Date: Jan 2010
Posts: 110
Default Using SRMA output to call variants with mpileup

I am trying to call variants with SAMTools mpileup using a BAM output of the SRMA program but am encountering problems whereby mpileup gives zero output.

It works fine when I use a BAM file generated pre-realignment so it is quite possible something is missing from my workflow post-realignment.

The original realignment was performed by splitting by chromosome so the subsequent steps I have employed are:

1) samtools merge -h original.sam realigned.bam *bam (merge all bam files)
2) samtools sort -n realigned.bam realigned_nsort.bam (sort by query names)
3) samtools fixmate realigned_nsort.bam realigned_fixed.bam (fix mate pair info)
4) MarkDuplicates I= realigned_fixed.bam O=final.bam (mark duplicate reads)
5) samtools sort realigned_fixed.bam final.bam (sort by coordinates)
6) samtools mpileup -uf genome.index final.bam

Mpileup generates nothing at this point. Can anyone suggest what might be going wrong?
gavin.oliver is offline   Reply With Quote
Old 07-19-2011, 05:56 AM   #2
gavin.oliver
Senior Member
 
Location: uk

Join Date: Jan 2010
Posts: 110
Default

I have discovered what the problem is.

Whilst SAMTools fixmate inserts pair mapping information that was removed by SRMA, it appears to mark all reads as having an unmapped mate. The reads are thus seen as anomalies by mpileup and ignored.

This seems to be happening because the sort -n command is not producing a properly sorted file.

E.g. the first three query names in the output file are

chr1_gl000191_random(Strand-Offset105824--106431)
chr1_gl000191_random(Strand+Offset5--434)
chr1_gl000191_random(Strand-Offset106058--106427)

The reads are paired so each sequence name should appear twice together.

Any ideas what could be causing the sort to go wrong?

Last edited by gavin.oliver; 07-19-2011 at 06:44 AM.
gavin.oliver is offline   Reply With Quote
Old 12-14-2011, 06:27 AM   #3
jkersh
Member
 
Location: boulder, co

Join Date: Apr 2010
Posts: 11
Default

I have just run into the same problem, has anyone found a way to work around the sort by name issue in Samtools? Or is there another way to do this?
jkersh is offline   Reply With Quote
Old 04-17-2014, 09:20 AM   #4
madonjoe
Junior Member
 
Location: San Diego

Join Date: Dec 2013
Posts: 6
Default Same problem encountered here

I also ran into the same problem here. I cannot do variant calling based on the newly re-aligned bam file. Is there a workaround?
madonjoe 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:09 AM.


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