This is a test version of Biostars. For the public version, visit https://www.biostars.org.
samtools mpileup keeps making an empty file

I'm trying to call SNPs from a pooled DNA sample using the PoPoolation2 pipeline. I'm stuck on a step that requires me to make a mpileup file from my BAM files. According to multiple tutorials this should work;

samtools mpileup -B input1.bam input2.bam > out.mpileup

However, when I try it makes an empty file.

I've tried adding a reference, re-sorting my BAM files and I checked for errors using Picardtool's ValidateSamFile. Even weirder, I can generate a file using GATK's Pileup command. Any ideas why samtools mpileup doesn't work for me?

samtools pileup
• 6,458 views • link
updated by anonym101 anonym101 • written by James Reeve James Reeve

version of samtools ? any error message ? number of reads per bam ? average mapping quality of the reads ? etc...

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Samtools v1.9. No error messages. Average 470 million reads per bam. MAPQ = 48.56.

• 0 views • link
written by James Reeve James Reeve

i do have a follow up question. Is using the -A option a good idea? I have same problem as OP (i.e. empty mpileup), but can't find any useful info. of course with -A it seems to work but...how reliable will all downstream analysis be? I asked the same question here but wondering whether you have other ideas...

• 0 views • link
written by anonym101 anonym101

2 answers

Try to pass -f <reference fasta> as a parameter before the bam files

• 0 views • link
written by gab • gab

mpileup should produce an output without reference.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I tired that, to no avail.

• 0 views • link
written by James Reeve James Reeve

show me two alignments from a random BAM please.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Alignment 1:

E00434:123:HM2NVCCXY:8:1101:8237:62452  1185    HiC_chromosome_1        1       60      34S117M =       47      197     GATCGTTGCTATATTTGTTCATAGTTTGTTTATTTGTTTTTCCACTGCTTTGCCAAAATAAACGGTAAAAAGTGTTTCCAAAACTTTTTCCTGGATATAACGCACACGATTTTCCAAGCACAAAAACATGGCGTATAACGGCAGGAGAGAC AAFFFAJJJ<JFJJJJJJJFJ7JJJJFJJJJFJ-FJJJJJAFAJ-FF<JJJJJJA7FJJJJJJJJJAJJJAJJ7FJAJ<AFAJJJFFJJAFFJJJFJAJ-AFJJJJ7A---AJJFAAF<JFFJAFFFAA-<A7AJ-----AA)77))-7-7 MD:Z:108C8      PG:Z:MarkDuplicates     RG:Z:Index_12   NM:i:1  AS:i:112        XS:i:72

Alignment 2:

E00434:123:HM2NVCCXY:8:2224:23906:72860 77  *   0   0   *   TTTAGTCTTTTTTGTCTTAATAGGCTCCCCTCGATTTAATTCACTCCCCCCAATTAATGAAAACAACTCAGCAGTTATCAACCCCATCAAACGACTTGCATGAGGGAGTATTCTCGCTGGTCTCTTAATTACTTCAGATCGGGAGAGCACC AAF<A<FFJJJJJFFFFJ-FFJ<JJFFJJJ<<--<JJJJJJ<-A7<FJFJJ-A<J<F--FJJJ-FJ--7<-7<7-7FJ<<J-AJJ-7<FJA7-A7-AAFFFAFJJJ<A--7<AJF<JJ7A--<A-<-<FFA7AA-<FJAFAF)F-AAFFJ) PG:Z:MarkDuplicates RG:Z:Index_12   AS:i:0  XS:i:0
• 0 views • link
written by James Reeve James Reeve

alignment 2 is not mapped, alignment 1 has sam flag 1185 = read is PCR or optical duplicate . Is there any non-exotic sam record in your sam file?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Here's one that is more normal

E00434:176:HV2G5CCXY:3:1109:27123:67955 97      chrI    1       60      50S79M  =       1       79      CTTCTCCTCTTCCTCTTCTTCCTCTTCCTCCTCCTCTTCTTCCTCTTCCTCTTCTTCTCCTCTTCCTCTTCCTTCTCCCTCTTCCTCCCTGTCGGCGTGTCATCAGATCTGACCAGTGTGTGTGTGTGT       AAFFFJJJJJJJJJJFJJJJJJJJJFJJJJJFJJJJJFJFJFJJJJJJJFFAJJJJFJJJJJFFJJJJJJJJJAF<FJJAJAAJFJFJJJJFJFJJ<7F<JFFJJJJJJJF<JJJFFJJ<7A7<J<FJF       SA:Z:chrXII,8525877,+,17S40M72S,0,2;    MD:Z:79 NM:i:0  AS:i:79 XS:i:27 RG:Z:IDT_64     PG:Z:MarkDuplicates
• 0 views • link
written by James Reeve James Reeve

this read is paired but "not properly paired" , pileup will ignore this

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I checked all my reads. Most reads are flagged as improperly paired. I toggled the option to not skip anomalous read pairs -A and it works. Thanks Pierre for walking me through this.

• 0 views • link
written by James Reeve James Reeve

Log in to answer this question.

More posts like this