samtools mpileup generates an empty output

[david.tamborero]

Hello,

I am runnning the samtools mpileup to a bam file generated by bfast (and processed by Picard RemoveDuplicates), and everything seems to work fine except for a certain sample in which the mpileup command generates an empty file.

The inputted bam file seems OK; it is an exome paired end alignement (of note both end pairs of each read are aligned, since I've found in the forum that it could be a cause of mpileup failure).

I've tried with several samtools versions and with several argument values, but it always rsults in the empty file after 1/2 hour of processing time.

I'm pretty confused with that, since I've used the same commands for all the samples and in general it worked fine. Anyone has a clue of what is happening?

Thanks!

[david.tamborero]

Just for those than can be interested in this thread, I've found that all the reads of this sample have been accounted as anomalous reads by mpileup, so they are excluded (and thus, an empty output is generated) if the -A flag of the mpileup command is not set.

cheers!
david

[amaer]

Thank you very much, I just ran into the exact same problem and I was completely puzzled.
I wonder why the mpileup considered all those reads anomalous, do you know what criteria it uses to make that call? Is it a problem with aligner used?

[david.tamborero]

I was checking this long time ago, and now I don't remember it on detail. I guess mpileup filtered out those reads that could be untrustwirthy, as those with zero-mapping quality or those improperly paired. For exact details, i guess you can ask to the samtools mailing list.

[ksowmya]

similar issue with RNA-Seq reads

Hi I have been trying to use samtools mpileup on RNA-Seq reads. These reads were in fastq format, and aligned against a mm10 Ensemble68 transcriptome reference using Bowtie2. After aligning I sorted the sam files by genome coordinates and then converted to bam. Now I am trying to use the same bam files and the reference to create the mpileup for further analysis but I keep getting an empty file.

mpileup only runs for about 2 secs and then stops. The only output I get is: [mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000

Does anyone have any suggestions?

[Richard Finney]

what is the samtools command you are using?

what does the header to the bam file you are using look like?

did you create an index for the bam file?

Does the output from samtools flagstat look good/

[ksowmya]

Hi Richard,

The command I am using is
samtools mpileup -f ./reference.fa ./D6-5066.bam > ./D6-5066-mpileup

When I created the sam file by aligning to the Ensemble68 reference, the sam file had a header, but I needed to use the command 'convert-iso-to-genome-coords' to convert the isoforms genome coordinates. The resulting sam does not have a header and so the bam file needed to be created using this, as recommended for sam files without @SQ header

samtools faidx ref.fa
samtools view -bt ref.fa.fai aln.sam > aln.bam

(link to info about command: http://dna.engr.uconn.edu/software/NGSTools/README.txt)

So the bam also does not have a header. Do you think these steps might be part of the problem? And no, I did not create an index for the bam files. There didn't seem to be any indication to do so, but I can see if it makes a difference.

I have not used flagstat before, and am not sure what it does.

[cwarden45]

I encountered this problem, which was caused by me accidentally using the forward read for both the R1 and R2 reads.

So, using "-A" might have allowed the .pileup file to be created without being empty or generating errors. However, I think it was important to understand why so many orphan "anomalous" reads were being flagged (which meant that I really needed to re-run the alignment).