Using SRMA output to call variants with mpileup

[gavin.oliver]

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]

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?

[jkersh]

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?

[madonjoe]

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?